set.seed(123)

1 Preliminary function definitions

In this section, we define some functions that we will use throughout the rest of this notebook. Each definition is hidden by default.

get_size_factor

get_size_factor <- function(cells) {
    coords <- cells@images$slice1@coordinates
    width <- max(coords$col) - min(coords$col)
    height <- max(coords$row) - min(coords$row)
    area <- width * height
    size_factor <- 140 / sqrt(area)

    return(size_factor)
}

rotate_image

A function to rotate the ST image by any angle.

rotate_image <- function(p, rot_angle) {
    if (is.null(rot_angle)) {
        rot_angle <- 0
    }
    gt <- ggplot_gtable(ggplot_build(p))
    panel_idx <- which(gt$layout$name == "panel")
    rot_vp <- viewport(angle = 360 - rot_angle)
    gt[["grobs"]][[panel_idx]] <- editGrob(gt[["grobs"]][[panel_idx]],
                                           vp = rot_vp)
    p_rot <- ggdraw() + draw_grob(gt)

    return(p_rot)
}

get_sample_name_from_plot

# Work out the sample name based on the number of spots in the plot.
get_sample_name_from_plot <- function(p) {
    # prepend "x" so that don't start with number
    sample_name <- switch(paste0("x", nrow(p$data)),
           x915 = "PCW_7_REP_3",
           x2523 = "PCW_11_REP_2",
           # (only 1855 spots if used for maturation path:)
           x1855 = "PCW_11_REP_2",
           # (only 2696 spots if used for maturation path:)
           x2696 = "PCW_17_REP_3",
           x3753 = "PCW_17_REP_3",
           x55 = "Donor1_day40_dif2_org2",
           x145 = "Donor2_day70_dif2_org2",
           x249 = "Donor1_day120_dif1_org1")

    return(sample_name)
}

FixSpatialPlot

FixSpatialPlot <- function(p, sample_name = NULL, num_cells = NULL) {
    # Run after all other modifications to Spatial*Plot (e.g. changes to theme,
    # scale colors, etc). Rotate the image and scale the dots appropriately. I
    # can't see a neater way of doing this for now (rebuilding the ggplot image
    # in `rotate_image` means that themes etc can no longer be applied:
    # stackoverflow.com/q/78027896/6914552).

    if (is.null(sample_name)) {
        sample_name <- get_sample_name_from_plot(p)
    }

    point_size_factors <- list("PCW_7_REP_3" = 2.4,
                               "PCW_11_REP_2" = 1.5,
                               "PCW_11_REP_2_sub" = 15,
                               "PCW_17_REP_3" = 1.8,
                               "Donor1_day40_dif2_org2" = 9,
                               "Donor2_day70_dif2_org2" = 7.5,
                               "Donor1_day120_dif1_org1" = 5.8,
                               "Donor1_day40_dif2_org2_sub" = 30,
                               "Donor2_day70_dif2_org2_sub" = 30,
                               "Donor1_day120_dif1_org1_sub" = 30)

    if ((!is.null(sample_name)) && (sample_name %in% names(point_size_factors))) {
        size_factor <- point_size_factors[[sample_name]]
    } else {
        size_factor <- 80 / sqrt(num_cells)
    }

    p$layers[[1]]$aes_params$point.size.factor <- size_factor
    p <- rotate_image(p, cells[[gsub("_sub", "", sample_name)]]@misc$rot_angle)

    return(p)
}

SpatialFeaturePlotBlend

gen_color_grid <- function(side_length, bottom_left, bottom_right, top_left,
                           top_right) {

    grad_gen <- function(start, end, n = side_length) {
        colfunc <- colorRampPalette(c(start, end))
        return(colfunc(n))
    }

    # x_y = "x to y"; "bl" = "bottom left", etc
    bl_tl <- grad_gen(bottom_left, top_left)
    br_tr <- grad_gen(bottom_right, top_right)

    l <- lapply(seq_len(length(bl_tl)),
           function(i) {
               start <- bl_tl[i]
               end <- br_tr[i]
               grad_gen(start, end)
           })

    return(matrix(unlist(l), ncol = side_length, nrow = side_length))
}

SpatialFeaturePlotBlend <- function(cells_obj, column_1, column_2,
                                    combine = TRUE, column_1_alt_name = NULL,
                                    column_2_alt_name = NULL,
                                    bottom_left = "#000000",
                                    bottom_right = "#0000FF",
                                    top_left = "#FF0000",
                                    top_right = "#FFFF00", ...)  {
    # Convert decimal number to hexadecimal. Pad with 0s if only a single
    # character following conversion.
    as_hex <- function(num) {
        hex_str <- as.character(as.hexmode(num))
        if (nchar(hex_str) == 1) {
            hex_str <- paste0("0", hex_str)
        }

        return(hex_str)
    }

    blend_plot_theme <- theme(legend.position = "none",
                              plot.title = element_text(hjust = 0.5))

    plot_list <- lapply(c(column_1, column_2),
                        function(column) {
                            max_color <- ifelse(column == column_1,
                                                "#FF0000", "#0000FF")
                            SpatialFeaturePlot(cells_obj, column, ...) +
                                scale_fill_gradient(low = "#000000",
                                                    high = max_color) +
                                ggtitle(column) +
                                blend_plot_theme
                        })

    dat <- FetchData(cells_obj, c(column_1, column_2)) %>%
            as.matrix()
    colnames(dat) <- c(column_1, column_2)
    side_length <- 100
    col_grid <- gen_color_grid(side_length, bottom_left, bottom_right,
                               top_left, top_right)
    dat_norm <- apply(dat, 2,
                      function(x) {
                          if (max(x) != 0) {
                              round((side_length - 1) * x / max(x)) + 1
                          } else {
                              rep(1, length(x))
                          }
                      })
    cols <- sapply(seq_len(nrow(dat_norm)),
                   function(x) {
                       col_grid[dat_norm[x, 1], dat_norm[x, 2]]
                   })

    new_md_column <- paste0(column_1, "_vs_", column_2)
    cells_obj[[new_md_column]] <- cols
    names(cols) <- as.character(cols)

    plot_list[[3]] <- SpatialDimPlot(cells_obj, new_md_column, cols = cols,
                                     ...) +
                        ggtitle(paste0(column_1, "_", column_2)) +
                        blend_plot_theme

    legend_grid <- expand.grid(seq(from = min(dat[, column_1]),
                                   to = max(dat[, column_1]),
                                   length.out = side_length),
                               seq(from = min(dat[, column_2]),
                                   to = max(dat[, column_2]),
                                   length.out = side_length))
    colnames(legend_grid) <- c(column_1, column_2)
    legend_colors <- c(col_grid)

    legend_grid$color <- legend_colors
    names(legend_colors) <- legend_colors

    legend <- ggplot(legend_grid,
                     aes(x = .data[[column_1]], y = .data[[column_2]],
                         color = color)) +
                geom_point(shape = 15, size = 1.9) +
                scale_color_manual(values = legend_colors) +
                coord_cartesian(expand = FALSE) +
                theme(legend.position = "none", aspect.ratio = 1,
                      panel.background = element_blank(),
                      axis.text.x = element_text(angle = 45, hjust = 1)) +
                xlab(ifelse(is.null(column_1_alt_name), column_1, column_1_alt_name)) +
                ylab(ifelse(is.null(column_2_alt_name), column_2, column_2_alt_name))

    plot_list[[4]] <- wrap_plots(ggplot() + theme_void(), legend,
                                 ggplot() + theme_void(), ncol = 1,
                                 heights = c(0.2, 0.6, 0.2))

    if (combine == FALSE) {
        return(plot_list)
    } else {
        p <- wrap_plots(plot_list, nrow = 1,
                        widths = c(0.28, 0.28, 0.28, 0.16))
        return(p)
    }
}

add_visual_reduction

add_visual_reduction <- function(cell_obj) {
    old_default_assay <- DefaultAssay(cell_obj)
    DefaultAssay(cell_obj) <- "Spatial"
    visual_coords <- GetTissueCoordinates(cell_obj)
    colnames(visual_coords) <- c("visual_1", "visual_2")
    vis_coords_mat <- as.matrix(visual_coords)
    cell_obj@reductions$visual <- CreateDimReducObject(key = "visual_",
                                                       embeddings = vis_coords_mat,
                                                       assay = "Spatial")
    DefaultAssay(cell_obj) <- old_default_assay

    return(cell_obj)
}

2 Load libraries

Load the R libraries required in the analysis. Also set some options for the generation of this document.

library(Seurat)
library(sctransform)
library(patchwork)
library(viridis)
library(ggplot2)
library(cowplot)
library(grid)
library(ggpubr)
library(ggrepel)
library(openxlsx)
library(DESeq2)
library(scuttle)
library(spacexr)
library(fgsea)
library(dplyr)
library(tibble)
library(reshape2)
library(ggtext)
library(knitr)
library(msigdbr)
library(showtext)
library(ggforce)

opts_chunk$set(message = FALSE, warning = FALSE, fig.align = "center",
               out.width = "75%")

Now set up some miscellaneous things that are important throughout:

# Set to true if you want to run even slow sections
compile_all <- TRUE

fig_list <- list("fig_5" = list(), "fig_6" = list(), "fig_7" = list())

main_figs_path <- "../paper_figs/main/"
supp_figs_path <- "../paper_figs/supplementary/"

fig_paths <- list("main" = main_figs_path, "supplementary" = supp_figs_path)

output_figure_format <- ".pdf"
# Width of A4 paper in cm
a4_width <- 21

# Override ggsave's new default black background color
ggsave <- function(plot_obj, ..., bg = "white") {
    print(plot_obj)
    ggplot2::ggsave(..., bg = bg)
}

# Set correct font family and minimum text size:
# First, load Arial using `showtext` library:
font_add("Arial", "arial_unicode.ttf")
font_add("ArialBold", "arial_bold.ttf")
showtext_auto() # This line is important!!
# Then override default theme and functions to enforce 11pt as a minimum and
# Arial as the font family:
# theme_grey is the default ggplot2 theme, and within it, the smallest text
# size appears to be 0.8 * base_size. So, by setting base_size to 13.75, this
# minimum becomes 11.
pub_fig_theme <- theme_grey(base_size = 13.75) +
                    theme(text = element_text(family = "Arial"),
                          panel.background = element_rect(fill = "white", colour = NA))
# and automatically apply this theme by default:
theme_set(pub_fig_theme)

# Overwrite ggplot2's annotate and geom_text functions:
annotate <- function(..., size = 11, family = "Arial") {
    ggplot2::annotate(..., size = size / .pt, family = family)
}
geom_text <- function(..., size = 11, family = "Arial") {
    ggplot2::geom_text(..., size = size / .pt, family = family)
}

# Abbreviated sample names
s_names_abbr <- list("PCW_7_REP_3" = "7 PCW", "PCW_11_REP_2" = "11 PCW",
                     "PCW_17_REP_3" = "17 PCW",
                     "Donor1_day40_dif2_org2" = "d40",
                     "Donor2_day70_dif2_org2" = "d70",
                     "Donor1_day120_dif1_org1" = "d120")

3 Load data

Load the datasets from file.

cells <- list()
samples <- c("OCT_BLOCK_1", "OCT_BLOCK_2", "OCT_BLOCK_3", "OCT_BLOCK_4",
             "PC_AGE_1_REP_1", "PC_AGE_1_REP_2", "PC_AGE_1_REP_3",
             "PC_AGE_2_REP_1", "PC_AGE_2_REP_2", "PC_AGE_2_REP_3",
             "PC_AGE_3_REP_1", "PC_AGE_3_REP_2", "PC_AGE_3_REP_3",
             "PC_AGE_4_REP_1", "PC_AGE_4_REP_2", "PC_AGE_4_REP_3")

for (s in samples) {
    dir_name <- paste0("../data/new_count_dirs/HGKFGDMXY_count_dirs/",
                       s, "/outs")
    # Load cell2location results for later use
    file_name <- paste0("../data/cell2loc_results_docker_3/cell2location_results_",
                        s, ".h5ad")
    h5ad <- zellkonverter::readH5AD(file_name)

    # Add timepoints to sample names
    if (grepl("PC_AGE", s)) {
        s <- sub("PC_AGE_1", "PCW_7", s)
        s <- sub("PC_AGE_2", "PCW_11", s)
        s <- sub("PC_AGE_3", "PCW_12", s)
        s <- sub("PC_AGE_4", "PCW_17", s)
    }

    cells[[s]] <- Load10X_Spatial(data.dir = dir_name) %>%
                  SCTransform(assay = "Spatial", verbose = FALSE) %>%
                  RunPCA() %>%
                  FindNeighbors(reduction = "pca")

    for (ct in h5ad@metadata$mod$factor_names) {
        ct_idx <- which(h5ad@metadata$mod$factor_names == ct)
        # Use 5th percentile of posterior distributions, since (as per the
        # documentation) this represents the value of cell abundance that the
        # model has high confidence in (aka ‘at least this amount is present’)
        # https://cell2location.readthedocs.io/en/latest/notebooks/cell2location_tutorial.html#Visualising-cell-abundance-in-spatial-coordinates

        cells[[s]][[paste0("c2l_", ct)]] <- h5ad@metadata$mod$post_sample_q05$w_sf[, ct_idx]
    }

    cells[[s]] <- add_visual_reduction(cells[[s]])

    message(paste("Loaded", s))
}

# Rotate main tissues
cells$PCW_7_REP_3@misc$rot_angle <- 180
cells$PCW_11_REP_2@misc$rot_angle <- 180
cells$PCW_17_REP_3@misc$rot_angle <- 270

# Rotate supplementary tissues
cells$PCW_7_REP_1@misc$rot_angle <- 167
cells$PCW_7_REP_2@misc$rot_angle <- 90
cells$PCW_11_REP_1@misc$rot_angle <- 142
cells$PCW_11_REP_3@misc$rot_angle <- 223
cells$PCW_17_REP_1@misc$rot_angle <- 264
cells$PCW_17_REP_2@misc$rot_angle <- 266

Later on, we will use lists of sample names that we want to include in the main or supplementary texts. We define these lists here.

nice_tissues <- c("PCW_7_REP_3", "PCW_11_REP_2", "PCW_17_REP_3")
nice_organoids <- c("Donor1_day40_dif2_org2", "Donor2_day70_dif2_org2",
                    "Donor1_day120_dif1_org1")
nice_samples <- c(nice_tissues, nice_organoids)

all_organoids <- names(cells)[grepl("Donor", names(cells))]

all_tissues_list <- list()
for (t in c(7, 11, 12, 17)) {
    pcw_name <- paste0("PCW_", t)
    all_tissues_list[[pcw_name]] <- paste0(pcw_name, paste0("_REP_", 1:3))
}

nice_tissues_list <- list("PCW_7" = "PCW_7_REP_3", "PCW_11" = "PCW_11_REP_2",
                          "PCW_17" = "PCW_17_REP_3")

all_organoids_list <- list()
for (donor_day in list(c(1, 40), c(1, 70), c(1, 120),
                       c(2, 40), c(2, 70), c(2, 120))) {
    as_str <- paste0("Donor", donor_day[1], "_day", donor_day[2])
    all_organoids_list[[as_str]] <- paste0(as_str,
                                           c("_dif1_org1", "_dif1_org2",
                                             "_dif2_org1", "_dif2_org2"))
}

supp_samples <- setdiff(as.character(c(unlist(all_tissues_list),
                                       unlist(all_organoids_list))),
                        nice_samples)

samples_list <- list("main" = nice_samples, "supplementary" = supp_samples)

4 Subset individial organoids

We are now in a position to extract the different organoids from each capture area. Since there are six organoids in each capture area, we must label them with the correct metadata. We do so by clustering the spatial co-ordinates of the spots into six clusters and then assigning the corresponding metadata.

m <- matrix(, nrow = 24, ncol = 3)
m[1, ] <- c("OCT_BLOCK_1", "0", "Donor2_day70_dif1_org2")
m[2, ] <- c("OCT_BLOCK_1", "1", "Donor1_day40_dif2_org1")
m[3, ] <- c("OCT_BLOCK_1", "2", "Donor2_day40_dif2_org2")
m[4, ] <- c("OCT_BLOCK_1", "3", "Donor2_day120_dif1_org2")
m[5, ] <- c("OCT_BLOCK_1", "4", "Donor2_day40_dif1_org1")
m[6, ] <- c("OCT_BLOCK_1", "5", "Donor1_day120_dif2_org2")
m[7, ] <- c("OCT_BLOCK_2", "0", "Donor2_day120_dif1_org1")
m[8, ] <- c("OCT_BLOCK_2", "1", "Donor1_day120_dif1_org2")
m[9, ] <- c("OCT_BLOCK_2", "2", "Donor2_day120_dif2_org2")
m[10, ] <- c("OCT_BLOCK_2", "3", "Donor1_day70_dif2_org2")
m[11, ] <- c("OCT_BLOCK_2", "4", "Donor1_day120_dif2_org1")
m[12, ] <- c("OCT_BLOCK_2", "5", "Donor1_day40_dif2_org2")
m[13, ] <- c("OCT_BLOCK_3", "0", "Donor1_day120_dif1_org1")
m[14, ] <- c("OCT_BLOCK_3", "1", "Donor2_day70_dif2_org2")
m[15, ] <- c("OCT_BLOCK_3", "2", "Donor2_day40_dif2_org1")
m[16, ] <- c("OCT_BLOCK_3", "3", "Donor2_day40_dif1_org2")
m[17, ] <- c("OCT_BLOCK_3", "4", "Donor2_day70_dif1_org1")
m[18, ] <- c("OCT_BLOCK_3", "5", "Donor1_day40_dif1_org1")
m[19, ] <- c("OCT_BLOCK_4", "0", "Donor1_day70_dif1_org1")
m[20, ] <- c("OCT_BLOCK_4", "1", "Donor1_day70_dif1_org2")
m[21, ] <- c("OCT_BLOCK_4", "2", "Donor2_day120_dif2_org1")
m[22, ] <- c("OCT_BLOCK_4", "3", "Donor2_day70_dif2_org1")
m[23, ] <- c("OCT_BLOCK_4", "4", "Donor1_day70_dif2_org1")
m[24, ] <- c("OCT_BLOCK_4", "5", "Donor1_day40_dif1_org2")
colnames(m) <- c("slide", "region", "sample_name")
df <- data.frame(m)

r <- data.frame(orig = c("dif1_org1", "dif1_org2", "dif2_org1", "dif2_org2"),
                replacement = paste0("n", as.character(1:4)))

for (s in c("OCT_BLOCK_1", "OCT_BLOCK_2", "OCT_BLOCK_3", "OCT_BLOCK_4")) {
    DefaultAssay(cells[[s]]) <- "Spatial"
    cells[[s]] <- add_visual_reduction(cells[[s]])

    cells[[s]] <- cells[[s]] %>%
                  FindNeighbors(dims = 1:2, reduction = "visual") %>%
                  FindClusters(resolution = 0.3, verbose = FALSE)
    if (s == "OCT_BLOCK_3") {
        # Combine clusters 4 and 6 for OCT_BLOCK_3 as they are the same organoid
        clust_6_idxs <- which(cells[[s]]$seurat_clusters == "6")
        cells[[s]]$seurat_clusters[clust_6_idxs] <- "4"
        cells[[s]]$seurat_clusters <- droplevels(cells[[s]]$seurat_clusters)
        Idents(cells[[s]]) <- cells[[s]]$seurat_clusters
    }

    cells[[s]]$org_name <- "unknown"
    for (i in levels(cells[[s]]$seurat_clusters)) {
        organoid_name <- subset(df, slide == s & region == i)$sample_name
        cells[[s]]$org_name[which(cells[[s]]$seurat_clusters == i)] <- organoid_name
        cells[[organoid_name]] <- subset(cells[[s]], seurat_clusters == i)
    }

    df_vis <- as.data.frame(cells[[s]]@reductions$visual@cell.embeddings)
    df_vis$seurat_clusters <- cells[[s]]$seurat_clusters
    new_visual_1 <- -1 * df_vis$visual_1
    # Flip y axis:
    df_vis$visual_1 <- new_visual_1 + (min(df_vis$visual_1) - min(new_visual_1))
    df_vis$ident <- NA

    cluster_means <- lapply(unique(df_vis$seurat_clusters),
                            function(x) {
                                subset(df_vis, seurat_clusters == x)[, -3] %>%
                                    apply(2, mean)
                            }) %>%
                        do.call(rbind, .)
    rownames(cluster_means) <- unique(df_vis$seurat_clusters)
    cluster_means <- as.data.frame(cluster_means)
    cluster_means$cluster_name <- rownames(cluster_means)
    cluster_means$ident <- NA
    for (i in cluster_means$cluster_name) {
        old_name <- subset(df, slide == s & region == i)$sample_name
        prefix <- substr(old_name, 1, nchar(old_name) - 9)
        old_suffix <- substr(old_name, nchar(old_name) - 8, nchar(old_name))
        new_suffix <- r[r$orig == old_suffix, "replacement"]
        new_name <- paste0(prefix, new_suffix)
        cluster_means[cluster_means$cluster_name == i, ]$cluster_name <- new_name
    }

    cluster_means_all_dots <- df_vis[, c("visual_1", "visual_2", "ident")]
    cluster_means_all_dots$cluster_name <- ""
    combined_df <- rbind(cluster_means, cluster_means_all_dots)
    p1 <- SpatialDimPlot(cells[[s]], pt.size.factor = 1, stroke = 0,
                         alpha = 0.6) +
            geom_point(data = df_vis, aes(x = visual_2, y = visual_1),
                       alpha = 0) +
            coord_cartesian(clip = "off") +
            geom_label_repel(data = combined_df,
                             aes(x = visual_2, y = visual_1,
                                 label = cluster_name),
                             max.overlaps = 10000, force = 50,
                             xlim = c(-Inf, Inf), ylim = c(-Inf, Inf),
                             fill = "white", direction = "both") +
            scale_x_continuous(expand = expansion(mult = 0.5)) +
            scale_y_continuous(expand = expansion(mult = 0.5)) +
            theme(legend.position = "none")

    p2 <- SpatialDimPlot(cells[[s]], alpha = 0) +
            theme(legend.position = "none")

    p <- ggarrange(p2, p1, nrow = 1, widths = c(1, 2.5))
    p <- annotate_figure(p, top = text_grob(s))

    fig_name <- paste0("organoid_names_", s)
    ggsave(p, filename = paste0(supp_figs_path, fig_name, output_figure_format),
           dpi = 300, width = 20, height = 10, units = "cm")

    print(p)
}

5 Remove outlier spots

for (s in c("PCW_7_REP_3", "PCW_11_REP_2", "Donor2_day70_dif2_org2")) {
    x <- cells[[s]]@reductions$visual@cell.embeddings[, 2]
    y <- cells[[s]]@reductions$visual@cell.embeddings[, 1]

    if (s == "PCW_7_REP_3") {
        is_outlier <- ((y + x) < 380) & (y < 180)
    } else if (s == "PCW_11_REP_2") {
        is_outlier <- (y < 150)
    } else if (s == "Donor2_day70_dif2_org2") {
        is_outlier <- ((0.8 * x) + y) > 578
    }

    cells[[s]][["is_outlier"]] <- is_outlier
    cells[[s]] <- subset(cells[[s]], is_outlier == FALSE)
}

6 Process main samples

# for (s in nice_samples) {
# organoids which cause a crash in the chunk below:
bad_organoids <- c("Donor2_day40_dif1_org1", "Donor1_day120_dif2_org2")
for (s in setdiff(names(cells), bad_organoids)) {
    message(s)
    cells[[s]] <- SCTransform(cells[[s]], assay = "Spatial",
                              verbose = FALSE) %>%
                    RunPCA() %>%
                    FindNeighbors(reduction = "pca", dims = 1:10)
}

7 Plot Cell2Location outputs

nice_ct_names <- list()
nice_ct_names[["c2l_Astro"]] <- "Astrocyte"
nice_ct_names[["c2l_Eryth"]] <- "Erythrocyte"
nice_ct_names[["c2l_OPC_1"]] <- "Oligodendrocyte precursor 1"
nice_ct_names[["c2l_OPC_2"]] <- "Oligodendrocyte precursor 2"
nice_ct_names[["c2l_Unk"]] <- "Unknown"
nice_ct_names[["c2l_VascLepto"]] <- "Vascular Lepto???"
nice_ct_names[["c2l_hDA1a"]] <- "Dopaminergic Neurons 1a"
nice_ct_names[["c2l_hDA1b"]] <- "Dopaminergic Neurons 1b"
nice_ct_names[["c2l_hDA2"]] <- "Dopaminergic Neurons 2"
nice_ct_names[["c2l_hDA3/hGABA/hSer"]] <- "Dopaminergic Neurons 3"
nice_ct_names[["c2l_hEndo"]] <- "Endothelial"
nice_ct_names[["c2l_hMgl"]] <- "Microglia"
nice_ct_names[["c2l_hMidPre"]] <- "Midbrain precursor"
nice_ct_names[["c2l_hNPro"]] <- "Neuron progenitor"
nice_ct_names[["c2l_hNbDA"]] <- "Dopaminergic neuroblast"
nice_ct_names[["c2l_hNbGaba"]] <- "GABA-ergic neuroblast"
nice_ct_names[["c2l_hPeric"]] <- "Pericyte"
nice_ct_names[["c2l_hPreDA"]] <- "Dopaminergic Neuron Precursor"
nice_ct_names[["c2l_hProgFPM"]] <- "Post-mitotic progentior???"
nice_ct_names[["c2l_hProgM"]] <- "Mitotic progenitor???"
nice_ct_names[["c2l_hRgl1"]] <- "Radial glial 1"
nice_ct_names[["c2l_hRgl2/immAstro"]] <- "Radial glial 2 / astrocyte???"
nice_ct_names[["c2l_hRgl3_caudal"]] <- "Radial glial 3 / caudal"
nice_ct_names[["c2l_hRgl4/MultiEpend"]] <- "Radial glial 4 / Ependyma???"

beautify_ct_plot <- function(p, plot_title = NULL) {
    p <- p + theme(legend.position = "right",
              legend.title = element_blank(),
              legend.text = element_blank(),
              plot.title = element_text(size = 8, hjust = 0.5)) +
        scale_fill_viridis_c() +
        guides(fill = guide_colorbar(ticks.colour = NA,
                                     barwidth = 0.5,
                                     barheight = 3)) +
        ggtitle(plot_title)
    p_spots <- p + theme(legend.position = "none")
    p_dat <- p$data
    c2l_col <- grep("c2l_", colnames(p_dat), value = TRUE)
    c2l_dat <- p_dat[, c2l_col]
    ct_max <- sprintf("%.1f", max(c2l_dat))
    p_leg <- as_ggplot(get_legend(p)) +
                    annotate(geom = "text", x = 0.45, y = 0.75,
                             label = ct_max) +
                    theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
    p_spots <- FixSpatialPlot(p_spots)
    p <- ggarrange(ggplot() + theme_void(), p_spots, p_leg,
                   ggplot() + theme_void(), nrow = 1,
                   widths = c(0.1, 0.6, 0.2, 0.1))

    return(p)
}
# Rename the slice's image to the sample name
for (s in nice_samples) {
    cells[[s]]@images[[s]] <- cells[[s]]@images[["slice1"]]
    cells[[s]]@images[["slice1"]] <- NULL
}

cells$merged_nice_samples <- merge(x = cells[[nice_samples[1]]],
                                   y = cells[nice_samples[2:length(nice_samples)]])

sample_names <- names(cells$merged_nice_samples@images)
num_samples <- length(sample_names)

custom_marker_plot <- function(m_plot, m_name) {
    m_plot <- m_plot +
                ggtitle(m_name) +
                theme(legend.position = "none",
                      plot.title = element_text(size = 7, hjust = 0.5,
                                                margin = margin(b = 0, unit = "cm"))) +
                guides(fill = guide_colorbar(label.vjust = -6,
                                             label.hjust = 1,
                                             label.position = "top",
                                             title.vjust = 0.2,
                                             ticks.colour = NA)) +
                scale_fill_gradient(low = "white", high = "red")

    return(m_plot)
}

marks_and_cts_list <- list()
marks_and_cts_list[[1]] <- list("c2l_hRgl1" = c("HES5", "OTX2"),
                                "c2l_hDA2" = c("TH", "KCNJ6"))
marks_and_cts_list[[2]] <- list("c2l_hNbDA" = c("DCX", "EBF2"),
                                "c2l_hDA1a" = c("TH", "LMO3"),
                                "c2l_hDA1b" = c("TH", "LMO3"),
                                "c2l_hDA2" = c("TH", "DDC"),
                                "c2l_hDA3/hGABA/hSer" = c("ALDH1A1", "KCNJ6"))
names(marks_and_cts_list) <- c("main", "supplementary")

for (ct_class in names(marks_and_cts_list)) {
    marks_and_cts <- marks_and_cts_list[[ct_class]]
    plot_list <- list()
    sample_names_plots <- lapply(s_names_abbr[Images(cells$merged_nice_samples)],
                                 function(sample_name) {
                                     ggplot() +
                                         annotate(geom = "text", x = 0, y = 0,
                                                  label = sample_name,
                                                  family = "ArialBold") +
                                         theme_void()
                                 })
    sample_names_plots_combined <- wrap_plots(sample_names_plots, nrow = 1)
    ct_name_width <- 0.8
    plot_list_widths <- c(ct_name_width, rep(1, num_samples))

    plot_list[["first_row"]] <- ggplot() + theme_void()
    plot_list[["sample_names"]] <- wrap_plots(ggplot() + theme_void(),
                                              sample_names_plots_combined,
                                              nrow = 1,
                                              widths = c(plot_list_widths[1],
                                                         sum(plot_list_widths[-1])))
    bottom_left <- "#000000"
    bottom_right <- "#0000FF"
    top_left <- "#FF0000"
    top_right <- "#FFFF00"

    for (ct in unique(names(marks_and_cts))) {
        to_plot <- c(ct, marks_and_cts[which(names(marks_and_cts) == ct)][[1]])
        # plot_list_inner arranged as: 1-6: CT plots; 7-12: M1 plots; 13-18: M2 plots
        to_plot_fix <- gsub("/", ".", to_plot)
        plot_list_inner <- SpatialFeaturePlot(cells$merged_nice_samples,
                                              to_plot_fix, stroke = NA,
                                              combine = FALSE)
        plot_list_inner_names <- apply(expand.grid(Images(cells$merged_nice_samples),
                                                   to_plot_fix),
                                       1,
                                       function(x) {
                                           paste0(x[1], "_", x[2])
                                       })
        names(plot_list_inner) <- plot_list_inner_names

        plot_list_mid <- list()
        ct_sub <- substr(ct, 5, nchar(ct))
        # Add newlines after dashes
        ct_sub <- gsub("/", "/\n", ct_sub)
        line_x <- 0.9
        plot_list_mid[["ct_name"]] <- ggplot() +
                                        xlim(-1, 1) +
                                        ylim(-1, 1) +
                                        annotate(geom = "text", x = 0, y = 0.5,
                                                 label = ct_sub, hjust = 0.5,
                                                 family = "ArialBold") +
                                        annotate(geom = "text", x = 0,
                                                 y = -0.4, label = to_plot[2],
                                                 hjust = 0.5, color = "red",
                                                 size = 8) +
                                        annotate(geom = "text", x = 0,
                                                 y = -0.6, label = to_plot[3],
                                                 hjust = 0.5, color = "blue",
                                                 size = 8) +
                                        geom_segment(aes(x = line_x,
                                                         xend = line_x,
                                                         y = -1, yend = 1)) +
                                        geom_segment(aes(x = line_x - 0.1,
                                                         xend = line_x,
                                                         y = 0.5,
                                                         yend = 0.5)) +
                                        theme_void()

        for (sample_idx in seq(num_samples)) {
            sample_name <- nice_samples[sample_idx]
            ct_plot <- plot_list_inner[[paste0(sample_name, "_",
                                               to_plot_fix[1])]]

            ct_max <- round(max(cells[[sample_name]][[to_plot[1]]]), 1)
            plot_margin <- unit(c(0, 0, 0, 0), "cm")
            ct_plot <- ct_plot +
                  theme(legend.position = "none",
                        legend.background = element_blank(),
                        legend.title = element_text(size = 5, hjust = 0),
                        legend.text = element_blank(),
                        plot.title = element_blank(),
                        plot.margin = plot_margin,
                        plot.background = element_blank()) +
                  scale_fill_viridis_c() +
                  labs(fill = paste(ct_max, "   ")) +
                  guides(fill = guide_colorbar(ticks.colour = NA,
                                               barwidth = 0.3,
                                               barheight = 1))

            ct_plot <- FixSpatialPlot(ct_plot)

            sfpbs <- SpatialFeaturePlotBlend(cells[[sample_name]],
                                             column_1 = to_plot[2],
                                             column_2 = to_plot[3],
                                             stroke = NA, combine = FALSE)
            marker_plots <- sfpbs[[3]] +
                                theme(plot.title = element_blank(),
                                      plot.margin = plot_margin,
                                      plot.background = element_blank())
            marker_plots <- FixSpatialPlot(marker_plots,
                                           sample_name = sample_name)

            p <- wrap_plots(ct_plot, marker_plots, ncol = 1)
            plot_list_mid[[sample_idx + 1]] <- p
        }
        plot_list[[ct]] <- wrap_plots(plot_list_mid, nrow = 1,
                                      widths = c(ct_name_width,
                                                 rep(1, num_samples)))
    }

    plot_list[["last_row"]] <- ggplot() + theme_void()

    heights_list <- list("main" = c(-0.4, 0.5, 4, 4, -1.3),
                         "supplementary" = c(-0.4, 0.5, 4, 4, 4, 4, 4, -1.3))
    p <- wrap_plots(plot_list, ncol = 1, heights = heights_list[[ct_class]])

    sfpb_leg <- sfpbs[[4]][[2]]
    x_min_max <- as.numeric(quantile(sfpb_leg$data[, 1], c(0, 1)))
    y_min_max <- as.numeric(quantile(sfpb_leg$data[, 2], c(0, 1)))

    sfpb_leg <- sfpb_leg +
                    theme(axis.ticks = element_blank(),
                          axis.title.x = element_text(color = bottom_right,
                                                      size = 8, vjust = 0.5),
                          axis.title.y = element_text(color = top_left,
                                                      size = 8, vjust = 0.5),
                          axis.text = element_text(size = 8)) +
                    xlab("Gene 1") +
                    ylab("Gene 2") +
                    scale_x_continuous(breaks = x_min_max,
                                       labels = c("Low", "High")) +
                    scale_y_continuous(breaks = y_min_max,
                                       labels = c("Low", "High"))

    p <- wrap_plots(p, sfpb_leg, nrow = 1, widths = c(0.9, 0.1))

    file_name <- paste0(fig_paths[[ct_class]], "c2l_combined_",
                        ct_class, output_figure_format)
    ggsave(p, filename = file_name, units = "cm", dpi = 300, width = a4_width,
           height = ifelse(ct_class == "main", 12, 30))

    if (ct_class == "main") {
        fig_list$fig_5[["c2l_combined_main"]] <- p
        fig_list$fig_5$row_1 <- p
    }
}

main_c2l <- c("c2l_hRgl1", "c2l_hProgFPM", "c2l_hNbDA", "c2l_hPreDA",
              "c2l_hDA1a", "c2l_hDA1b", "c2l_hDA2", "c2l_hDA3/hGABA/hSer")
supp_c2l_1 <- c("c2l_hRgl2/immAstro", "c2l_hRgl3_caudal",
                "c2l_hRgl4/MultiEpend", "c2l_hMidPre", "c2l_hNPro",
                "c2l_hProgM", "c2l_hNbGaba", "c2l_Astro")
supp_c2l_2 <- c("c2l_OPC_1", "c2l_OPC_2", "c2l_hMgl", "c2l_hEndo",
                "c2l_hPeric", "c2l_Eryth", "c2l_VascLepto", "c2l_Unk")
c2l_cts <- list("main" = list("main" = main_c2l),
                "supplementary" = list("supplementary_1" = supp_c2l_1,
                                       "supplementary_2" = supp_c2l_2))

strip_prefix <- function(x) {
    substr(x, 5, nchar(x))
}

c2l_cts_pretty <- lapply(c2l_cts, function(x) {lapply(x, strip_prefix)})
names(c2l_cts_pretty$main$main) <- main_c2l
names(c2l_cts_pretty$supplementary$supplementary_1) <- supp_c2l_1
names(c2l_cts_pretty$supplementary$supplementary_2) <- supp_c2l_2

fig_paths <- list("main" = main_figs_path, "supplementary" = supp_figs_path)

sample_group <- "supplementary"
for (ct_group in names(c2l_cts[[sample_group]])) {
    cts <- c2l_cts[[sample_group]][[ct_group]]
    plot_list_outer <- list()
    pretty_ct_names <- c2l_cts_pretty[[sample_group]][[ct_group]]
    ct_names_plots <- lapply(cts,
                             function(ct) {
                                 ggplot() +
                                     annotate("text", x = 0, y = 0, size = 8,
                                              label = pretty_ct_names[[ct]],
                                              angle = 90) +
                                     theme_void()
                             })
    plot_list_outer[["ct_names"]] <- wrap_plots(ct_names_plots, ncol = 1)
    num_cts <- length(cts)
    num_samples <- length(nice_samples)
    for (s in nice_samples) {
        plot_list <- lapply(cts,
                            function(cell_type) {
                                p <- SpatialFeaturePlot(cells[[s]],
                                                        cell_type,
                                                        image.alpha = 0) +
                                        theme(legend.position = "none",
                                              plot.title = element_blank(),
                                              plot.margin = unit(c(0, 0, 0, 0), "cm")) +
                                    scale_fill_viridis_c()
                                p <- FixSpatialPlot(p, sample_name = s, num_cells = ncol(cells[[s]]))
                                if (cell_type == cts[1]) {
                                    p_s <- ggplot() +
                                     annotate("text", x = 0, y = 0,
                                              size = 8,
                                              label = s_names_abbr[[s]]) +
                                     theme_void()
                                    p <- wrap_plots(p_s, p, ncol = 1,
                                                    heights = c(0.1, 0.9))
                                }
                                return(p)
                            })
        plot_list_outer[[s]] <- wrap_plots(plot_list, ncol = 1,
                                           nrow = num_cts)
    }
    title_percentage <- 5
    widths <- c(title_percentage,
                 rep((100 - title_percentage) / num_samples, num_samples))
    p <- wrap_plots(plot_list_outer, nrow = 1, ncol = num_samples + 1,
                    widths = widths)
    f_name <- paste0(fig_paths[[sample_group]], "c2l_combined_",
                     ct_group, "_no_markers", output_figure_format)
    ggsave(p, filename = f_name, dpi = 300, width = a4_width,
           height = 29.7, units = "cm", limitsize = FALSE)
}

p <- SpatialFeaturePlot(cells[[s]], c2l_cts[[sample_group]][[ct_group]][1]) +
        scale_fill_viridis_c(breaks = c(0.05, 0.33),
                             labels = c("Low", "High")) +
        theme(legend.title = element_blank()) +
        guides(fill = guide_colourbar(ticks.colour = NA))
p <- cowplot::get_legend(p) %>% as_ggplot()
ggsave(p,
       filename = paste0(fig_paths[["supplementary"]],
                         "c2l_supp_legend", output_figure_format),
       dpi = 300, units = "cm", height = 2, width = 4)

8 RNA molecule plots

braun_genes <- c("STMN2", "NKX2-1", "NKX2-2", "BARHL1", "SHH", "EN2", "EMX2",
                 "WNT5A")
colors_to_use <- c("#FFBFCC", "#77AADD", "#EE8866", "#F2DC6B", "#99DDFF",
                   "#44BB99", "#BBCC33", "#8C8C8C")
names(colors_to_use) <- braun_genes

generate_molecule_plot <- function(obj, genes, cols) {
    ad <- GetAssayData(obj, assay = "Spatial", layer = "counts")
    p <- SpatialFeaturePlot(obj, features = "TH", alpha = 0)
    p_data <- p[[1]]$data

    df <- data.frame(matrix(ncol = 0, nrow = 0))
    for (g in genes) {
        gt_0_spots <- names(which(ad[g, ] > 0))
        gt_0_expr <- ad[g, gt_0_spots]
        gt_0_spots_duped <- rep(gt_0_spots, gt_0_expr)
        if (length(gt_0_spots) > 0) {
            new_rows <- cbind(p_data[gt_0_spots_duped, c("imagecol", "imagerow")],
                              gene = g, "TH" = 0)
            df <- rbind(df, new_rows)
        }
    }

    df$imagerow_flip <- max(df$imagerow) + (-1 * df$imagerow)
    df$gene <- factor(df$gene, levels = genes)
    p <- SpatialFeaturePlot(obj, features = "TH", alpha = 0,
                            image.alpha = 0.5) +
    geom_point(data = df, aes(x = imagecol, y = imagerow_flip,
                              color = gene),
               position = position_jitter(width = 3, height = 3,
                                          seed = 123),
               alpha = 1, size = 0, stroke = 0.2, show.legend = TRUE) +
    scale_color_manual(values = cols, drop = FALSE) +
    guides(fill = "none",
           color = guide_legend(override.aes = list(size = 3), ncol = 2)) +
    theme(legend.title = element_blank(),
          legend.position = "right",
          legend.key = element_rect(fill = NA))

    return(p)
}

p_for_legend <- generate_molecule_plot(cells[[nice_samples[1]]], braun_genes,
                                       colors_to_use)
leg <- as_ggplot(get_legend(p_for_legend))

arrow_plot <- ggplot() +
                geom_segment(aes(x = 0, xend = 0, y = -0.25, yend = 0.25),
                             arrow = arrow(length = unit(0.25, "cm"),
                                           ends = "both",
                                           type = "closed")) +
                annotate("text", x = 0, y = 0.4, label = "D") +
                annotate("text", x = 0, y = -0.4, label = "V") +
                scale_x_continuous(limits = c(-1, 1)) +
                scale_y_continuous(limits = c(-1, 1)) +
                theme_void()

plot_list <- list()
for (s in nice_tissues) {
    p <- generate_molecule_plot(cells[[s]], braun_genes, colors_to_use)
    p <- p + theme(legend.position = "none")
    p <- rotate_image(p, cells[[s]]@misc$rot_angle) %>%
            annotate_figure(top = text_grob(s_names_abbr[[s]],
                                            size = 11,
                                            family = "ArialBold",
                                            vjust = 1.5,
                                            hjust = 0.3))
    plot_list[[s]] <- p
}
p <- ggarrange(plotlist = plot_list, nrow = 1)
p_complete <- ggarrange(ggplot() + theme_void(), arrow_plot, p, leg,
                        ggplot() + theme_void(), nrow = 1,
                        widths = c(3, 1, 14, 7, 3))

fig_list$fig_5$molecules <- p_complete
fig_list$fig_5$row_4 <- p_complete

f_name <- paste0(fig_paths[["main"]], "molecule_plots", output_figure_format)
ggsave(p_complete, filename = f_name, width = a4_width - 3,
       height = 4, units = "cm")

9 Supergroups

Carry out clustering and then plot the expression some general marker genes with respect to these clusters

nice_samples <- nice_samples[grep(pattern = "OCT_", x = nice_samples, invert = TRUE)]
marks_list <- list("Progenitor" = c("HES1", "SOX2", "VIM"),
                   "Post-mitotic" = c("TH", "DCX", "MAP2"))
num_clusters <- 2

cluster_res_list <- list()
cluster_res_list[["PCW_7_REP_3"]] <- list("2" = 0.05, "3" = 0.1)
cluster_res_list[["PCW_11_REP_2"]] <- list("2" = 0.045, "3" = 0.08)
cluster_res_list[["PCW_17_REP_3"]] <- list("2" = 0.04, "3" = 0.05)
cluster_res_list[["Donor1_day40_dif2_org2"]] <- list("2" = 0.3)
cluster_res_list[["Donor2_day70_dif2_org2"]] <- list("2" = 0.4)
cluster_res_list[["Donor1_day120_dif1_org1"]] <- list("2" = 0.2)

plot_list_spatial <- list()
plot_list_dotplot <- list()
for (s in nice_samples) {
    DefaultAssay(cells[[s]]) <- "SCT"
    res <- cluster_res_list[[s]][[as.character(num_clusters)]]
    if (!is.null(res)) {
        cells[[s]] <- FindClusters(cells[[s]], resolution = res)
    } else {
        res <- 0.01
        cells[[s]] <- FindClusters(cells[[s]], resolution = res)
        while (length(levels(cells[[s]]$seurat_clusters)) < 2) {
            res <- res + 0.01
            cells[[s]] <- FindClusters(cells[[s]], resolution = res)
        }
    }

    cells[[s]]$cluster_ident <- as.character(cells[[s]]$seurat_clusters)
    zero_idxs <- as.numeric(which(cells[[s]]$cluster_ident == "0"))
    avg_vim_zeroes <- mean(GetAssayData(cells[[s]])["VIM", zero_idxs])
    avg_vim_non_zeroes <- mean(GetAssayData(cells[[s]])["VIM", -zero_idxs])

    if ((!is.nan(avg_vim_zeroes)) && (avg_vim_zeroes > avg_vim_non_zeroes)) {
        cells[[s]]$cluster_ident[zero_idxs] <- "Progenitor"
        cells[[s]]$cluster_ident[-zero_idxs] <- "Post-mitotic"
    } else {
        cells[[s]]$cluster_ident[zero_idxs] <- "Post-mitotic"
        cells[[s]]$cluster_ident[-zero_idxs] <- "Progenitor"
    }
    cells[[s]]$cluster_ident <- factor(cells[[s]]$cluster_ident,
                                       levels = names(marks_list))
    Idents(cells[[s]]) <- cells[[s]]$cluster_ident

    cluster_color_list <- list("Progenitor" = "#5B9EDA",
                               "Post-mitotic" = "#D05280")

    marks_color_list <- lapply(names(marks_list),
                               function(x) {
                                   rep(cluster_color_list[[x]],
                                       length(marks_list[[x]]))
                               })
    marks_colors <- unlist(marks_color_list)

    p1 <- SpatialDimPlot(cells[[s]], stroke = NA) +
            scale_fill_manual(values = lapply(levels(Idents(cells[[s]])),
                                              function(x) {
                                                  cluster_color_list[[x]]
                                              })) +
            theme(legend.position = "none")
    p1 <- FixSpatialPlot(p1, sample_name = s, num_cells = ncol(cells[[s]]))
    plot_list_spatial[[s]] <- p1 %>% annotate_figure(top = text_grob(s_names_abbr[[s]],
                                                                     size = 11,
                                                                     family = "ArialBold",
                                                                     vjust = 1.5))

    y_axis_order <- levels(Idents(cells[[s]]))
    cluster_colors <- as.character(unlist(cluster_color_list[y_axis_order]))
        cols_to_use <- c("#FFFFFF", "#000000")
    p2 <- DotPlot(cells[[s]], features = as.character(unlist(marks_list)),
                  assay = "Spatial", scale = FALSE, dot.scale = 4) +
            theme(axis.title = element_blank(),
                  legend.title = element_blank(),
                  strip.text = element_blank(),
                  axis.text.x = element_markdown(angle = 90, hjust = 1,
                                                 color = marks_colors, size = 8),
                  axis.text.y = element_markdown(size = 11,
                                                 color = cluster_colors),
                  legend.key.height = unit(0.2, 'cm'),
                  plot.margin = unit(c(0, 0, 0, 0), "cm")) +
            ylab("Cluster") +
            scale_y_discrete(limits = y_axis_order)
    max_val <- max(p2$data$avg.exp.scaled)
    min_val <- min(p2$data$avg.exp.scaled)
    p2 <- p2 + scale_color_gradientn(colours = c("#FFFFFF", "#000000"),
                                     breaks = (max_val - min_val) * c(0.1, 0.9),
                                     labels = c("Low", "High")) +
            guides(color = guide_colourbar(ticks.colour = NA, barwidth = 0.4))

    if (s != nice_samples[1]) {
        p2 <- p2 + theme(axis.text.y = element_blank(),
                         axis.line.y = element_blank(),
                         axis.ticks.y = element_blank())
    }

    if (s != nice_samples[length(nice_samples)]) {
        p2 <- p2 + theme(legend.position = "none")
    }

    plot_list_dotplot[[s]] <- p2
}
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 915
## Number of edges: 25770
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9539
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2523
## Number of edges: 82152
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9574
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 3753
## Number of edges: 121879
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9633
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 55
## Number of edges: 1124
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7101
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 145
## Number of edges: 5931
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6270
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 249
## Number of edges: 8060
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8116
## Number of communities: 2
## Elapsed time: 0 seconds
spatial_widths <- rep(1, length(plot_list_spatial))
spatial_plots <- wrap_plots(plot_list_spatial, nrow = 1, widths = spatial_widths)

dotplot_widths <- rep(1, length(plot_list_dotplot))
dotplot_widths[1] <- 1
dotplot_widths[length(dotplot_widths)] <- 1
dotplots <- wrap_plots(plot_list_dotplot, nrow = 1, widths = dotplot_widths) & theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))

p <- wrap_plots(spatial_plots, dotplots, ncol = 1, heights = c(0.8, 0.3))

fig_name <- "supergroups"

fig_list$fig_5$supergroups <- p
fig_list$fig_5$row_3 <- p
ggsave(p, filename = paste0(main_figs_path, fig_name,
                            output_figure_format),
       dpi = 300, width = a4_width, height = 8, units = "cm", limitsize = FALSE)

10 Differential Expression with CSIDE

The CSIDE software allows us to identify differentially expressed genes with respect to some covariate (in our case, maturation). We run it below, first loading the cell2location deconvolution results.

10.1 Create reference object

pau_cells <- readRDS("../data/mlo_resolution075_Annot.RDS")

cluster_names <- c(0:23)
names(cluster_names) <- c("hRgl2/immAstro", "hNbDA", "hProgFPM", "OPC_1",
                          "VascLepto", "hDA1b", "hRgl1", "hDA1a",
                          "hRgl3_caudal", "hDA2", "hProgM", "hPreDA",
                          "hMidPre", "hMgl", "hEndo", "hNbGaba", "hNPro",
                          "hDA3/hGABA/hSer", "Unk", "hRgl4/MultiEpend",
                          "Astro", "hPeric", "Eryth", "OPC_2")
pau_cells$seurat_clusters_24_Annot <- names(cluster_names[pau_cells$seurat_clusters])

# Following https://satijalab.org/seurat/articles/spatial_vignette.html#spatial-deconvolution-using-rctd
ref <- pau_cells
ref <- UpdateSeuratObject(ref)
counts <- ref[["RNA"]]$counts
ref$seurat_clusters_24_Annot <- gsub("/", "_", ref$seurat_clusters_24_Annot)
Idents(ref) <- "seurat_clusters_24_Annot"
cluster <- as.factor(ref$seurat_clusters_24_Annot)

names(cluster) <- colnames(ref)
nUMI <- ref$nCount_RNA
names(nUMI) <- colnames(ref)
reference <- Reference(counts, cluster, nUMI)

rm(pau_cells)
rm(ref)
rm(counts)
gc()
##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells  14230070  760.0   21900860  1169.7   21900860  1169.7
## Vcells 849270702 6479.5 4178214862 31877.3 5175077528 39482.8

Now we are ready to run CSIDE.

convert_to_spatial_rna <- function(s) {
    coords <- GetTissueCoordinates(s)
    colnames(coords) <- c("x", "y")
    coords[is.na(colnames(coords))] <- NULL
    counts <- s[["Spatial"]]$counts

    return(SpatialRNA(coords, counts, colSums(counts)))
}

run_cside <- function(cell_obj, reference, sample_name, exp_var_name = "maturation_path", exp_var_id_1 = NULL) {
    RCTD <- create.RCTD(convert_to_spatial_rna(cell_obj),
                        reference = reference, max_cores = 1,
                        CELL_MIN_INSTANCE = 15)

    # The two lines are slightly hacky ways around bugs in CSIDE. Cell2Location
    # effectively runs in  full mode (in the terminology of RCTD)
    RCTD@config$RCTDmode <- "full"
    RCTD@config$doublet_mode <- FALSE

    c2l_weights <- cell_obj@meta.data[, grepl("^c2l", colnames(cell_obj@meta.data))]
    c2l_weights <- normalize_weights(c2l_weights)
    colnames(c2l_weights) <- gsub("/", "_", colnames(c2l_weights))
    colnames(c2l_weights) <- gsub("c2l_", "", colnames(c2l_weights))

    RCTD <- import_weights(RCTD, c2l_weights)
    RCTD <- set_cell_types_assigned(RCTD)

    if (exp_var_name == "maturation_path") {
        exp_var <- cell_obj$maturation
    } else if (exp_var_name == "mat_region") {
        exp_var <- as.integer(cell_obj[["mat_region"]] == exp_var_id_1)
    }
    names(exp_var) <- colnames(cell_obj)

    # Setting weight_threshold and cell_type_threshold much lower than their
    # default values in CSIDE, otherwise it crashes. This change is probably
    # required because cell2location is likely to generate much smaller weights
    # than RCTD as it assumes more cells per spot.
    cside_out <- run.CSIDE.single(RCTD, exp_var,
                                  doublet_mode = FALSE,
                                  weight_threshold = 0.2,
                                  cell_type_threshold = 1,
                                  fdr = 0.05)

    if (!is.null(sample_name)) {
        # Generate combined spreadsheet
        res <- cside_out@de_results$all_gene_list
        res <- res[which(unlist(lapply(res, nrow)) > 0)]
        for (ct in names(res)) {
            res[[ct]] <- cbind(gene = rownames(res[[ct]]), res[[ct]],
                               cell_type = ct)
            is_sig <- res[[ct]]$gene %in% rownames(cside_out@de_results$sig_gene_list[[ct]])
            res[[ct]]$is_signif <- is_sig
        }
        cside_out_merged <- do.call(rbind, res)
    }

    return(cside_out_merged)
}

We now define the maturaion covariates according to the anatomy of the tissue slices.

# Define PCW_11_REP_2 maturation

cells$PCW_11_REP_2 <- add_visual_reduction(cells$PCW_11_REP_2)

vis_emb <- cells$PCW_11_REP_2@reductions$visual@cell.embeddings
# Need to shift coords here to centre the axes at the correct point
x_coord <- vis_emb[, 2] - 255
y_coord <- max(vis_emb[, 1]) - vis_emb[, 1]
y_coord <- y_coord - 100
xy_old <- cbind(x_coord, y_coord)
cells$PCW_11_REP_2$x_coord <- x_coord
cells$PCW_11_REP_2$y_coord <- y_coord

# Need to do some rotating
rotation_angle <- 20 * (pi / 180) # 20 degrees counterclockwise (in radians)
rot_mat <- matrix(c(cos(rotation_angle), -sin(rotation_angle),
                    sin(rotation_angle), cos(rotation_angle)),
                  ncol = 2)
xy_new <- rot_mat %*% t(xy_old)
cells$PCW_11_REP_2$x_new <- xy_new[1, ]
cells$PCW_11_REP_2$y_new <- xy_new[2, ]

cells$PCW_11_REP_2$theta_new <- ifelse(cells$PCW_11_REP_2$y_new > 0,
                                       atan(abs(cells$PCW_11_REP_2$x_new / cells$PCW_11_REP_2$y_new)),
                                       pi - atan(abs(cells$PCW_11_REP_2$x_new / cells$PCW_11_REP_2$y_new)))

cells_sub <- subset(cells$PCW_11_REP_2, y_coord >= 0 & y_coord <= 280)

in_first_mat_traj <- (abs(cells_sub$x_new) <= 20) & (cells_sub$y_new >= -10)
prop_first <- prop.table(table(in_first_mat_traj))[["TRUE"]]

cells_sub$first_mat_traj <- cells_sub$y_new
cells_sub$second_mat_traj <- cells_sub$theta_new

norm_0_1 <- function(x) {
    x <- x - min(x)
    x <- x / max(x)

    return(x)
}

# Make two trajectories join and go from 0 -> 1 smoothly
first_mat_traj_renorm <- prop_first * norm_0_1(cells_sub$first_mat_traj[in_first_mat_traj])
second_mat_traj_renorm <- prop_first + ((1 - prop_first) * norm_0_1(cells_sub$second_mat_traj[!in_first_mat_traj]))

cells_sub$maturation <- NA
cells_sub$maturation[in_first_mat_traj] <- first_mat_traj_renorm
cells_sub$maturation[!in_first_mat_traj] <- second_mat_traj_renorm
cells$PCW_11_REP_2_sub <- cells_sub

if (compile_all) {
    cside_outs <- list()
    cside_outs[["PCW_11_REP_2"]] <- run_cside(cells_sub, reference, "PCW_11_REP_2")
}
## 
##            Astro            Eryth            hDA1a            hDA1b 
##              189               88             2843             3509 
##             hDA2  hDA3_hGABA_hSer            hEndo             hMgl 
##             2587              374              914             1242 
##          hMidPre            hNbDA          hNbGaba            hNPro 
##             1577             5340              689              542 
##           hPeric           hPreDA         hProgFPM           hProgM 
##              147             2021             4486             2269 
##            hRgl1   hRgl2_immAstro     hRgl3_caudal hRgl4_MultiEpend 
##             3359             6408             2688              215 
##            OPC_1            OPC_2              Unk        VascLepto 
##             3815               15              350             3617
# Define maturation for PCW_17_REP_3

cells$PCW_17_REP_3 <- add_visual_reduction(cells$PCW_17_REP_3)

vis_emb <- cells$PCW_17_REP_3@reductions$visual@cell.embeddings
x_coord <- vis_emb[, 2]
y_coord <- max(vis_emb[, 1]) - vis_emb[, 1]
y_coord <- y_coord
xy_old <- cbind(x_coord, y_coord)
cells$PCW_17_REP_3$x_coord <- x_coord
cells$PCW_17_REP_3$y_coord <- y_coord
cells$PCW_17_REP_3$y_coord_1 <- y_coord - 250
cells_sub <- subset(cells$PCW_17_REP_3, x_coord >= 145)

cells_sub$theta_1 <- atan((max(cells_sub$x_coord) - cells_sub$x_coord) / cells_sub$y_coord_1)
cells_sub$exclude <- (cells_sub$theta_1 <= 0) & (cells_sub$theta_1 >= -0.5)
cells_sub <- subset(cells_sub, exclude == FALSE)
cells_sub <- subset(cells_sub,
                    (y_coord < 100) & (x_coord < 250) & (c2l_hDA2 < 0.4),
                    invert = TRUE)

cells_sub$use_first_mat_traj <- cells_sub$y_coord >= 270
prop_first <- prop.table(table(cells_sub$use_first_mat_traj))[["TRUE"]]

first_mat_traj_idxs <- cells_sub$use_first_mat_traj
cells_sub$maturation <- NA
cells_sub$first_maturation <- norm_0_1(-cells_sub$x_coord)
cells_sub$maturation[first_mat_traj_idxs] <- prop_first * norm_0_1(cells_sub$first_maturation[first_mat_traj_idxs])

cells_sub$theta <- (pi / 2) - atan((max(cells_sub$x_coord) - cells_sub$x_coord) / (max(cells_sub$y_coord) - cells_sub$y_coord))
cells_sub$maturation[first_mat_traj_idxs == FALSE] <- prop_first + ((1 - prop_first) * norm_0_1(cells_sub$theta[first_mat_traj_idxs == FALSE]))
cells$PCW_17_REP_3_sub <- cells_sub

if (compile_all) {
    cside_outs[["PCW_17_REP_3"]] <- run_cside(cells_sub, reference, "PCW_17_REP_3")
}
## 
##            Astro            Eryth            hDA1a            hDA1b 
##              189               88             2843             3509 
##             hDA2  hDA3_hGABA_hSer            hEndo             hMgl 
##             2587              374              914             1242 
##          hMidPre            hNbDA          hNbGaba            hNPro 
##             1577             5340              689              542 
##           hPeric           hPreDA         hProgFPM           hProgM 
##              147             2021             4486             2269 
##            hRgl1   hRgl2_immAstro     hRgl3_caudal hRgl4_MultiEpend 
##             3359             6408             2688              215 
##            OPC_1            OPC_2              Unk        VascLepto 
##             3815               15              350             3617
# Define maturation for PCW_7_REP_3

cells$PCW_7_REP_3 <- add_visual_reduction(cells$PCW_7_REP_3)

vis_emb <- cells$PCW_7_REP_3@reductions$visual@cell.embeddings
x_coord <- vis_emb[, 2]
y_coord <- max(vis_emb[, 1]) - vis_emb[, 1]
xy_old <- cbind(x_coord, y_coord)
cells$PCW_7_REP_3$x_coord <- x_coord
cells$PCW_7_REP_3$y_coord <- y_coord
cells$PCW_7_REP_3$maturation <- norm_0_1(cells$PCW_7_REP_3$y_coord)

if (compile_all) {
    cside_outs[["PCW_7_REP_3"]] <- run_cside(cells$PCW_7_REP_3, reference, "PCW_7_REP_3")
}
## 
##            Astro            Eryth            hDA1a            hDA1b 
##              189               88             2843             3509 
##             hDA2  hDA3_hGABA_hSer            hEndo             hMgl 
##             2587              374              914             1242 
##          hMidPre            hNbDA          hNbGaba            hNPro 
##             1577             5340              689              542 
##           hPeric           hPreDA         hProgFPM           hProgM 
##              147             2021             4486             2269 
##            hRgl1   hRgl2_immAstro     hRgl3_caudal hRgl4_MultiEpend 
##             3359             6408             2688              215 
##            OPC_1            OPC_2              Unk        VascLepto 
##             3815               15              350             3617
# Define maturation for organoids
for (s in nice_organoids) {
    mins <- colMins(cells[[s]]@reductions$visual@cell.embeddings,
                    useNames = TRUE)
    maxs <- colMaxs(cells[[s]]@reductions$visual@cell.embeddings,
                    useNames = TRUE)
    center <- c(mins[1] + (maxs[1] - mins[1]) / 2,
                mins[2] + (maxs[2] - mins[2]) / 2)
    coords <- cells[[s]]@reductions$visual@cell.embeddings
    cells[[s]]$maturation <- apply(coords, 1, function(x) {
                                       sqrt((x[1] - center[1])^2 + (x[2] - center[2])^2)
                })
    cells[[s]]$maturation <- norm_0_1(cells[[s]]$maturation)

    if (compile_all) {
        cside_outs[[s]] <- run_cside(cells[[s]], reference, s)
    }
}
## 
##            Astro            Eryth            hDA1a            hDA1b 
##              189               88             2843             3509 
##             hDA2  hDA3_hGABA_hSer            hEndo             hMgl 
##             2587              374              914             1242 
##          hMidPre            hNbDA          hNbGaba            hNPro 
##             1577             5340              689              542 
##           hPeric           hPreDA         hProgFPM           hProgM 
##              147             2021             4486             2269 
##            hRgl1   hRgl2_immAstro     hRgl3_caudal hRgl4_MultiEpend 
##             3359             6408             2688              215 
##            OPC_1            OPC_2              Unk        VascLepto 
##             3815               15              350             3617 
## 
##            Astro            Eryth            hDA1a            hDA1b 
##              189               88             2843             3509 
##             hDA2  hDA3_hGABA_hSer            hEndo             hMgl 
##             2587              374              914             1242 
##          hMidPre            hNbDA          hNbGaba            hNPro 
##             1577             5340              689              542 
##           hPeric           hPreDA         hProgFPM           hProgM 
##              147             2021             4486             2269 
##            hRgl1   hRgl2_immAstro     hRgl3_caudal hRgl4_MultiEpend 
##             3359             6408             2688              215 
##            OPC_1            OPC_2              Unk        VascLepto 
##             3815               15              350             3617 
## 
##            Astro            Eryth            hDA1a            hDA1b 
##              189               88             2843             3509 
##             hDA2  hDA3_hGABA_hSer            hEndo             hMgl 
##             2587              374              914             1242 
##          hMidPre            hNbDA          hNbGaba            hNPro 
##             1577             5340              689              542 
##           hPeric           hPreDA         hProgFPM           hProgM 
##              147             2021             4486             2269 
##            hRgl1   hRgl2_immAstro     hRgl3_caudal hRgl4_MultiEpend 
##             3359             6408             2688              215 
##            OPC_1            OPC_2              Unk        VascLepto 
##             3815               15              350             3617
if (compile_all) {
    for (s in names(cside_outs)) {
        cside_outs[[s]]$sample_name <- s
    }
    out_spreadsheet <- createWorkbook()
    df <- do.call(rbind, cside_outs)
    sheet_name <- "maturation_degs"
    addWorksheet(out_spreadsheet, sheet_name)
    writeData(out_spreadsheet, sheet_name, df, rowNames = FALSE)
    saveWorkbook(out_spreadsheet, "../out_data/all_maturation_degs.xlsx",
                 overwrite = TRUE)
}
plot_list <- list()
for (sample_name in nice_samples) {
    if (sample_name == "PCW_11_REP_2") {
        s <- "PCW_11_REP_2_sub"
    } else if (sample_name == "PCW_17_REP_3") {
        s <- "PCW_17_REP_3_sub"
    } else {
        s <- sample_name
    }

    line_width <- 0.5
    end_arrow <- arrow(type = "closed", angle = 30,
                       length = unit(0.05, units = "npc"))
    p <- SpatialFeaturePlot(cells[[s]], "maturation", stroke = NA) +
            scale_fill_gradient(low = "yellow", high = "red",
                                breaks = c(0.1, 0.9),
                                labels = c("Low", "High"))

    if (sample_name == "PCW_7_REP_3") {
        p <- p + geom_segment(aes(x = p$data$imagecol[which.min(p$data$maturation)],
                                  y = p$data$imagerow[which.max(p$data$maturation)] + 10),
                              xend = p$data$imagecol[which.max(p$data$maturation)],
                              yend = p$data$imagerow[which.min(p$data$maturation)] - 10,
                              arrow = end_arrow,
                              linewidth = line_width)
    } else if (sample_name == "PCW_11_REP_2") {
        first_seg_end_x <- 200
        first_seg_end_y <- 310
        min_y <- min(p$data$imagerow)
        p <- p + geom_segment(aes(x = 250, y = 195.5,
                                  xend = first_seg_end_x, yend = first_seg_end_y),
                              linewidth = line_width) +
                 geom_curve(aes(x = first_seg_end_x, y = first_seg_end_y,
                                xend = 150, yend = 210),
                            linewidth = line_width,
                            curvature = 0.9, ncp = 10, angle = 30,
                            arrow = end_arrow) +
                 geom_curve(aes(x = first_seg_end_x, y = first_seg_end_y,
                                xend = 400, yend = 210),
                            linewidth = line_width,
                            arrow = end_arrow, curvature = -0.9, ncp = 10,
                            angle = 140)
    } else if (sample_name == "PCW_17_REP_3") {
        first_seg_end_x <- 200
        first_seg_end_y <- 380
        p <- p + geom_segment(aes(x = 450, y = 380, xend = first_seg_end_x,
                                  yend = first_seg_end_y),
                              linewidth = line_width) +
                 geom_curve(aes(x = first_seg_end_x, y = first_seg_end_y + 1,
                                xend = 380, yend = 200),
                            linewidth = line_width,
                            arrow = end_arrow, curvature = 0.5, ncp = 10,
                            angle = 90)
    } else {
        # Organoid
        # For some reason, finding the co-ordinates more programatically seems
        # to break things
        if (sample_name == "Donor1_day40_dif2_org2") {
            p <- p + geom_segment(aes(x = 333, y = 215, xend = 368,
                                      yend = 215),
                                  linewidth = line_width, arrow = end_arrow)
        } else if (sample_name == "Donor2_day70_dif2_org2") {
            p <- p + geom_segment(aes(x = 402, y = 188, xend = 450,
                                      yend = 188),
                                  linewidth = line_width, arrow = end_arrow)
        } else {
            p <- p + geom_segment(aes(x = 285, y = 309, xend = 330,
                                      yend = 309),
                                  linewidth = line_width, arrow = end_arrow)
        }
    }

    if (substr(sample_name, 1, 3) == "PCW") {
        plot_title <- paste(strsplit(sample_name, "_")[[1]][2], "pcw")
    } else {
        d <- strsplit(substr(sample_name, 11, nchar(sample_name)), "_")[[1]][1]
        plot_title <- paste0("d", d)
    }

    p_leg <- p + theme(legend.position = "left") +
                labs(fill = "Maturation")
    p <- p + theme(legend.position = "none",
                   plot.margin = unit(c(0, 0, 0, 0), "cm"))

    p <- FixSpatialPlot(p, sample_name = sample_name)
    p <- p + annotate(geom = "text", x = 0.5, y = 1, hjust = 0.5,
                      label = plot_title, size = 8)

    plot_list[[sample_name]] <- p
}
tissue_plots <- ggarrange(plotlist = plot_list[1:3], nrow = 1)
organoid_plots <- ggarrange(plotlist = plot_list[4:6], nrow = 1)
p <- ggarrange(ggplot() + theme_void(), tissue_plots, organoid_plots,
               ggplot() + theme_void(), nrow = 4, heights = c(0.15, 1, 1, 0))

fig_list$fig_5[["maturation_paths"]] <- p

file_name <- paste0(fig_paths[["main"]], "maturation_paths",
                    output_figure_format)
ggsave(p, filename = file_name, units = "cm", dpi = 300, height = 6,
       width = 10)

10.2 Plotting CSIDE results

window_length <- 50
max_idx <- max(sapply(nice_samples,
                      function(s) {
                          s <- ifelse(s %in% nice_tissues[2:3], paste0(s, "_sub"), s)
                          num_spots <- ncol(cells[[s]])
                          return(num_spots - window_length)
                      }))

plot_list <- list()
for (gene in c("TH", "CAMK2N1", "LMO4", "SPOCK1", "MTRNR2L12", "MT-ND5")) {
        mat_df <- data.frame(matrix(ncol = 3, nrow = 0))
        colnames(mat_df) <- c("Sample", "gene_expr", "mat_idx")
        for (sample_name in nice_samples) {
            s <- ifelse(sample_name %in% nice_tissues[2:3],
                        paste0(sample_name, "_sub"), sample_name)
            ad <- GetAssayData(cells[[s]], assay = "SCT")
            mat_order <- order(cells[[s]]$maturation)
            g_expr_wrt_mat <- ad[which(rownames(ad) == gene), ][mat_order]
            avg_expr <- sapply(1:(length(g_expr_wrt_mat) - window_length),
                               function(x) {
                                   mean(g_expr_wrt_mat[x:(x + window_length)])
                               })
            avg_expr <- avg_expr / max(avg_expr)
            num_spots_minus_window <- ncol(ad) - window_length
            mat_idx <- max_idx * seq(num_spots_minus_window) / num_spots_minus_window
            new_rows <- data.frame("Sample" = sample_name,
                                   "gene_expr" = avg_expr,
                                   "mat_idx" = mat_idx)
            mat_df <- rbind(mat_df, new_rows)
        }

    x_arrow <- arrow(type = "closed", angle = 30,
                     length = unit(0.5, units = "npc"))
    colors <- c("#FFC2C2", "#C36767", "#FF0000", "#9B9AE4", "#4592D4",
                "#000EFF")

    mat_df$Sample <- factor(mat_df$Sample, levels = nice_samples)
    num_lines <- length(unique(mat_df$Sample))
    expr_range <- max(mat_df$gene_expr) - min(mat_df$gene_expr)

    p <- ggplot(mat_df,
                aes(x = mat_idx, y = gene_expr, color = Sample)) +
            geom_line() +
            theme_minimal() +
            theme(panel.grid = element_blank(),
                  legend.title = element_blank(),
                  axis.title = element_text(size = 6),
                  axis.line.x = element_line(arrow = x_arrow),
                  axis.line.y = element_line(arrow = x_arrow),
                  axis.text = element_blank(),
                  plot.title = element_blank(),
                  legend.margin = margin(t = -0.3, l = 0.3, unit = "cm")) +
            xlab("Maturation") +
            ylab(gene) +
            scale_color_manual(values = colors[seq(1, num_lines)], labels = s_names_abbr) +
            guides(colour = guide_legend(ncol = 1, override.aes = list(linewidth = 1))) +
            scale_x_continuous(breaks = c(0.05, 0.95) * max(mat_df$mat_idx)) +
            scale_y_continuous(breaks = min(mat_df$gene_expr) + (c(0.05, 1) * expr_range)) +
            theme(plot.margin = unit(c(0, 0, 0.3, 0.3), "cm"))
            pub_fig_theme

     plot_list[[gene]] <- p
}

p <- ggarrange(plotlist = plot_list, nrow = 2, ncol = 3, common.legend = TRUE,
               legend = "right") +
        theme(legend.title = element_blank())
fig_list$fig_5[["all_genes_avg_per_sample"]] <- p
ggsave(p, filename = paste0(fig_paths[["supplementary"]],
                            "all_genes_avg_per_sample", output_figure_format),
       dpi = 300, width = 40, height = 12, units = "cm")

fig_list$fig_5$row_2 <- ggarrange(ggplot() + theme_void(),
                                  fig_list$fig_5$maturation_paths,
                                  fig_list$fig_5$all_genes_avg_per_sample,
                                  ggplot() + theme_void(),
                                  nrow = 1, widths = c(-0.5, 1, 3, -0.3)) &
                            pub_fig_theme &
                            theme(axis.title = element_blank())
fig_list$fig_5$row_2 <- ggarrange(ggplot() + theme_void(),
                                  fig_list$fig_5$row_2,
                                  ggplot() + theme_void(),
                                  ncol = 1, heights = c(0, 1, -0.1))

10.3 DEGs between maturation regions

flip <- function(x) {
    # Flip y axis of SpatialDimPlot
    min_x <- min(x)
    max_x <- max(x)
    at_0 <- x - min_x
    flipped <- -1 * at_0
    return(flipped + min_x + max_x)
}

markers <- c("TH", "CALB1", "ALDH1A1")
s <- "PCW_17_REP_3"
dot_plots <- list()
vln_plots <- list()
spatial_plots <- list()
expr_plots <- list()
cells_to_use <- cells[[s]]

region_colors <- list("1" = "#D81B60", "2" = "#1E88E5", "3" = "#FFC107",
                      "4" = "#78ECD9")
genes_to_plot <- c("EN1", "TH", "DDC", "TUBA1A", "NNAT", "SLC6A4", "FEV",
                   "SLC17A6", "NRN1", "ADCY1", "NR2F1", "SCG2")

ad <- GetAssayData(cells_to_use, slot = "data")[markers, ]
min_ad <- min(ad)
max_ad <- max(ad)

maturation_ids <- read.csv(paste0("../data/maturation_labellings/", s,
                                  "_mat_regions_TH.csv"))
maturation_ids_idxs <- match(colnames(cells_to_use),
                             maturation_ids$Barcode)
cells_to_use$mat_region <- maturation_ids[maturation_ids_idxs,
                                          "subcluster"]

non_th56_idxs <- (!cells_to_use$mat_region %in% c("TH_5", "TH_6")) | is.na(cells_to_use$mat_region)
cells_to_use <- cells_to_use[, non_th56_idxs]

cells_to_use$mat_region_num <- sapply(cells_to_use$mat_region,
                                      function(x) {
                                      substr(x, nchar(x), nchar(x) + 1)})

cells_sub <- cells_to_use[, !is.na(cells_to_use$mat_region_num)]

df <- cbind(GetTissueCoordinates(cells_sub), cells_sub$mat_region_num)
colnames(df) <- c("imagerow", "imagecol", "mat_region_num")
df[, markers] <- 1
df$imagerow <- flip(df$imagerow)

# Compute convex hulls around each group of points:
chulls <- data.frame(matrix(ncol = ncol(df), nrow = 0))
colnames(chulls) <- colnames(df)
for (i in unique(df$mat_region_num)) {
    df_sub <- subset(df, mat_region_num == i)
    chull_df_sub <- df_sub[chull(df_sub[, c("imagerow", "imagecol")]), ]
    chulls <- rbind(chulls, chull_df_sub)
}
breaks <- round(range(c(min_ad, max_ad)))

for (g in markers) {
    fp <- SpatialFeaturePlot(cells_to_use, g, stroke = NA) +
            scale_fill_gradient(low = "white", high = "red",
                                limits = c(min_ad, max_ad),
                                breaks = seq(breaks[1], breaks[2], 1)) +
            labs(fill = "Expr") +
            theme(legend.position = "none")

    p <- fp
    expansion_amount <- 0.006
    for (region in unique(chulls$mat_region_num)) {
        p <- p +
                geom_shape(data = subset(chulls, mat_region_num == region),
                           aes(x = imagecol,
                               y = imagerow + 110),
                           fill = NA,
                           color = region_colors[[as.character(region)]],
                           expand = expansion_amount)
    }

    p <- FixSpatialPlot(p, sample_name = s)

    p <- p %>%
            annotate_figure(top = text_grob(g))

    spatial_plots[[g]] <- p

    all_pairs <- t(combn(sort(unique(cells_to_use$mat_region_num)), 2))

    fm_outs_list <- list()
    for (i in seq(1, nrow(all_pairs))) {
        pair <- all_pairs[i, ]
        message(pair)
        marks <- FindMarkers(cells_to_use, group.by = "mat_region_num",
                             ident.1 = pair[1], ident.2 = pair[2],
                             verbose = FALSE)
        marks$gene <- rownames(marks)
        pair_name <- paste0(pair[1], "_", pair[2])
        marks$comparison <- pair_name
        fm_outs_list[[pair_name]] <- marks
    }

    mat_region_degs <- do.call(rbind, fm_outs_list) %>%
                        filter(p_val_adj <= 0.05) %>%
                        arrange(-avg_log2FC) %>%
                        distinct(gene) %>%
                        pull(gene)

    cells_sub <- cells_to_use[, is.na(cells_to_use$mat_region_num) == FALSE]

    dot_plots[[g]] <- DotPlot(cells_sub, features = mat_region_degs,
                              group.by = "mat_region_num", dot.scale = 4) +
                        theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
                              axis.title.x = element_blank(),
                              legend.text = element_text(size = 8),
                              legend.title = element_text(size = 8)) +
                        scale_size_continuous(breaks = c(0, 50, 100)) +
                        scale_color_viridis(breaks = c(-1, 0, 1)) +
                        scale_y_discrete(limits = rev) +
                        labs(color = "Avg. expression") +
                        ylab("TH region") +
                        guides(color = guide_colorbar(barwidth = 1,
                                                      barheight = 2))

    # Select only spots positive for g
    g_pos_idxs <- GetAssayData(cells_sub)[g, ] > 0
    cells_sub <- cells_sub[, g_pos_idxs]
    ad <- GetAssayData(cells_sub[genes_to_plot, ],
                       slot = "data")
    breaks <- floor(range(ad))
    df <- data.frame(t(ad), mat_region_num = cells_sub$mat_region_num,
                     check.names = FALSE)
    df <- tidyr::pivot_longer(df, cols = genes_to_plot,
                              values_to = "value", names_to = "gene")
    p <- ggplot(df, aes(x = mat_region_num, y = value,
                        fill = mat_region_num)) +
            geom_violin(linewidth = 0) +
            theme(legend.position = "none",
                  panel.grid.major.y = element_line(size = 0.1,
                                                    color = "#ABABAB"),
                  panel.grid.minor.y = element_line(size = 0.1,
                                                    color = "#ABABAB")) +
            geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf,
                          ymax = Inf),
                      colour = "black", fill = NA, inherit.aes = FALSE) +
            scale_y_continuous(breaks = breaks[2],
                               minor_breaks = seq(breaks[1], breaks[2],
                                                  1)) +
            scale_fill_manual(values = region_colors) +
            facet_wrap(vars(gene), ncol = 1,
                       strip.position = "right") +
            xlab(paste0(g, "+ spots")) +
            ylab("Expression") +
            theme(axis.line.x = element_line(size = 0.2),
                  axis.ticks.x = element_blank(),
                  axis.title = element_text(size = 10),
                  panel.spacing = unit(0, 'lines'),
                  strip.background = element_rect(fill = NA),
                  strip.text.y = element_text(angle = 0, size = 8))
    vln_plots[[g]] <- p
}

leg_df <- data.frame(x = 1:4, y = 1:4, mat_region_num = names(region_colors))
p_leg_l <- ggplot(leg_df, aes(x = x, y = y, color = mat_region_num)) +
            geom_point() +
            scale_color_manual(values = region_colors) +
            labs(color = "Region")
p_leg_l <- as_ggplot(get_legend(p_leg_l))

p_leg_r <- as_ggplot(get_legend(fp + theme(legend.position = "right")))

all_spatial_plots <- ggarrange(plotlist = spatial_plots, nrow = 1)
all_spatial_plots <- ggarrange(p_leg_l, all_spatial_plots, p_leg_r, nrow = 1,
                               widths = c(0.1, 0.8, 0.1))
fig_list$fig_6$row_1 <- all_spatial_plots

vln_plots <- ggarrange(plotlist = vln_plots, nrow = 1)
all_expr_plots <- ggarrange(dot_plots[["TH"]], vln_plots, ncol = 1,
                            heights = c(0.35, 0.65))
fig_list$fig_6$row_2 <- all_expr_plots

11 Cell-to-cell communication

We use the COMMOT software to evaluate the expression of ligand-receptor pairs in our tissue. The tool is implemented in Python, so this first chunk has already been run for the three samples of interest

samples_for_commot <- c("PCW_11_REP_2", nice_organoids)

for (s in samples_for_commot) {
    as_sce <- as.SingleCellExperiment(cells[[s]], assay = "Spatial")
    out_file <- paste0("../data/", s, "_for_commot.h5ad")
    zellkonverter::writeH5AD(sce = as_sce, file = out_file)
}
# Need to load this package so that we can use the COMMOT docker image in the
# chunk below.
library(docknitr)

# Notice that, in the options to the chunk below, we bodge the image name to
# allow us to mount the correct path to the data. This is necessary because we
# will run the chunk in a different docker to the one that the rest of the
# notebook is running in, and we therefore need to pass the absolute path to
# the data on the host machine (i.e. outside the docker container), rather than
# to the mounted location, which is what docknitr mounts by default.
import commot as ct
import scanpy as sc
import pandas as pd
import numpy as np

df_ligrec = ct.pp.ligand_receptor_database(database="CellChat", species="human")
sample_list = [
    "PCW_11_REP_2",
    "Donor1_day40_dif2_org2",
    "Donor2_day70_dif2_org2",
    "Donor1_day120_dif1_org1",
]

for sample_name in sample_list:
    print(sample_name)
    adata = sc.read_h5ad(f"/workdir2/data/{sample_name}_for_commot.h5ad")
    adata.obsm["spatial"] = adata.obsm["VISUAL"]
    adata.var_names_make_unique()
    sc.pp.normalize_total(adata, inplace=True)
    sc.pp.log1p(adata)

    # Only look at pairs of cells at max distance 1000um (i.e. 10 spots)
    ct.tl.spatial_communication(
        adata,
        database_name="user_database",
        df_ligrec=df_ligrec,
        dis_thr=1000,
        heteromeric=True,
    )

    for db in ["sender", "receiver"]:
        out_name = f"/workdir2/out_data/{sample_name}_commot_{db}_db.csv"
        adata.obsm[f"commot-user_database-sum-{db}"].to_csv(out_name)
## PCW_11_REP_2
## Donor1_day40_dif2_org2
## Donor2_day70_dif2_org2
## Donor1_day120_dif1_org1
nice_sr_pairs <- c("PTN.PTPRZ1", "MDK.PTPRZ1", "MDK.LRP1", "MIF.ACKR3")
gen_blended_plot <- function(cells_obj, sr_pair,
                             s_name, r_name,
                             bottom_right_col,
                             show_image = TRUE) {
    pl <- SpatialFeaturePlotBlend(cells_obj,
                                  paste0("s.", sr_pair),
                                  paste0("r.", sr_pair),
                                  combine = FALSE,
                                  column_1_alt_name = s_name,
                                  column_2_alt_name = r_name,
                                  bottom_right = bottom_right_col,
                                  stroke = 0,
                                  image.alpha = ifelse(show_image, 1, 0))
    pl[[3]] <- pl[[3]] + theme(plot.title = element_blank())

    return(pl[[3]])
}
bottom_right_col <- "#00FCFF"

# Assume square. In form: c(top_left_x, top_left_y, side_length)
rect_positions <- list("PCW_11_REP_2" = list(c(280, 200, 40),
                                             c(230, 260, 40),
                                             c(230, 310, 40)),
                       "Donor1_day40_dif2_org2" = list(c(310, 195, 20)),
                       "Donor2_day70_dif2_org2" = list(c(360, 160, 20)),
                       "Donor1_day120_dif1_org1" = list(c(300, 280, 20)))
sr_pair <- "PTN.PTPRZ1"

outer_plot_list <- list()
outer_plot_list_no_zooms <- list()
for (s in samples_for_commot) {
    file_prefix <- paste0("../out_data/", s, "_commot_")
    sender_db <- read.csv(paste0(file_prefix, "sender_db.csv"))
    receiver_db <- read.csv(paste0(file_prefix, "receiver_db.csv"))

    rownames(sender_db) <- sender_db[, 1]
    rownames(receiver_db) <- receiver_db[, 1]

    commot_matrix <- t(cbind(sender_db[, -1], receiver_db[, -1]))
    old_row_names <- rownames(commot_matrix)
    commot_matrix <- apply(commot_matrix, 2, as.numeric)
    rownames(commot_matrix) <- old_row_names
    # Remove "total.total" rows:
    commot_matrix <- commot_matrix[-which(grepl("total.total", rownames(commot_matrix))), ]
    commot_matrix <- commot_matrix[, match(colnames(cells[[s]]), colnames(commot_matrix))]

    commot_assay <- CreateAssayObject(data = commot_matrix, key = "COMMOT_")
    cells[[s]][["COMMOT"]] <- commot_assay

    cells[[s]] <- FindSpatiallyVariableFeatures(object = cells[[s]],
                                                assay = "COMMOT",
                                                slot = "data",
                                                selection.method = "moransi")
    spat_var <- SpatiallyVariableFeatures(cells[[s]], assay = "COMMOT",
                                          selection.method = "moransi")
    spat_var_without_prefix <- sapply(spat_var,
                                      function(x) {
                                          paste0(strsplit(x, ".", fixed = TRUE)[[1]][-1],
                                                 collapse = ".")
                                      })
    uniq_spat_var_pairs <- unique(spat_var_without_prefix)

    ad <- GetAssayData(cells[[s]], assay = "COMMOT")
    num_spots_expr <- function(sr) {
        sum(apply(ad[grepl(sr, rownames(ad)), ], 2,
                  function(x) {
                      any(x > 0.1)
                  }))
    }
    highly_expressed_idxs <- sapply(uniq_spat_var_pairs, num_spots_expr) > 10
    uniq_spat_var_pairs <- uniq_spat_var_pairs[highly_expressed_idxs]

    sr_pairs_to_plot <- head(c(nice_sr_pairs, uniq_spat_var_pairs), n = 9)

    old_default_assay <- DefaultAssay(cells[[s]])
    DefaultAssay(cells[[s]]) <- "COMMOT"

    all_s_0 <- all(ad[paste0("s.", sr_pair), ] == 0)
    all_r_0 <- all(ad[paste0("r.", sr_pair), ] == 0)
    if (all_s_0 | all_r_0) {
        return(ggplot() + theme_void())
    }
    splat <- strsplit(sr_pair, ".", fixed = TRUE)
    s_name <- splat[[1]][1]
    r_name <- splat[[1]][2]

    bp <- gen_blended_plot(cells[[s]], sr_pair, s_name, r_name,
                           bottom_right_col)
    bp_out <- bp

    border_colors <- c("#0000FF", "#CC33FF", "#CFEB99")

    bp_zoom_list <- list()
    for (rect_pos in rect_positions[[s]]) {
        rect_pos_str <- paste0(rect_pos, collapse = "_")
        rect_x_min <- rect_pos[1]
        rect_x_max <- rect_x_min + rect_pos[3]
        rect_y_min <- rect_pos[2]
        rect_y_max <- rect_y_min + rect_pos[3]
        border_color <- border_colors[length(bp_zoom_list) + 1]

        bp_out <- bp_out + geom_rect(xmin = rect_x_min, xmax = rect_x_max,
                                         ymin = rect_y_min, ymax = rect_y_max,
                                     color = border_color, fill = NA,
                                     linewidth = 0.75)

        bp_zoom_list[[rect_pos_str]] <- gen_blended_plot(cells[[s]], sr_pair,
                                                         s_name, r_name,
                                                         bottom_right_col,
                                                         show_image = FALSE) +
                                        coord_cartesian(xlim = c(rect_x_min,
                                                                 rect_x_max),
                                                        ylim = c(rect_y_min,
                                                                 rect_y_max)) +
                                        theme(panel.border = element_rect(colour = border_color,
                                                                          fill = NA,
                                                                          size = 2))
        bp_zoom_list[[rect_pos_str]] <- FixSpatialPlot(bp_zoom_list[[rect_pos_str]],
                                                   paste0(s, "_sub"))
    }

    p <- ggarrange(FixSpatialPlot(bp_out),
                   ggplot() + theme_void(),
                   ggarrange(plotlist = bp_zoom_list, ncol = 1),
                   nrow = 1,
                   widths = c(ifelse(length(bp_zoom_list) > 1, 1.5, 0.5),
                              ifelse(length(bp_zoom_list) > 1, 0, -0.1),
                              0.5))

    outer_plot_list[[s]] <- p
    ggsave(p,
           filename = paste0("../paper_figs/main/commot_nice_and_spat_var_sr_pairs_", s,
                             ".pdf"),
           dpi = 300, units = "cm", height = 20, width = 20)

    splat <- strsplit(sr_pair, ".", fixed = TRUE)
    s_name <- splat[[1]][1]
    r_name <- splat[[1]][2]
    pl <- lapply(sr_pairs_to_plot,
                 function(sr_pair) {
                     p <- SpatialFeaturePlotBlend(cells[[s]], paste0("s.", sr_pair),
                                             paste0("r.", sr_pair),
                                             combine = FALSE,
                                             bottom_right = bottom_right_col,
                                             stroke = 0)[[3]] +
                            ggtitle(gsub(".", " → ", sr_pair, fixed = TRUE)) +
                                        theme(plot.title = element_text(hjust = 0.5, size = 7))

                     FixSpatialPlot(p)
                 })
    outer_plot_list_no_zooms[[s]] <- ggarrange(plotlist = pl, nrow = 1)

    DefaultAssay(cells[[s]]) <- old_default_assay
}

p <- ggarrange(plotlist = outer_plot_list_no_zooms, ncol = 1)
ggsave(p, filename = paste0("../paper_figs/main/commot_nice_and_spat_var_sr_",
                            "pairs_no_zooms", output_figure_format),
       dpi = 300, units = "cm", height = 15, width = 30)

p <- ggarrange(outer_plot_list[[1]],
               ggarrange(plotlist = outer_plot_list[2:4], ncol = 1),
               ncol = 2, widths = c(0.6, 0.4))

p_leg <- SpatialFeaturePlotBlend(cells$PCW_7_REP_3, "nFeature_Spatial",
                                 "nCount_Spatial", combine = FALSE,
                                 column_1_alt_name = "Ligand",
                                 column_2_alt_name = "Receptor",
                                 bottom_right = bottom_right_col,
                                 stroke = 0)[[4]][[2]]
legend_break_percentiles <- c(0.1, 0.9)
p_leg <- p_leg +
            scale_x_continuous(breaks = as.numeric(quantile(p_leg$data[, 1],
                                                            legend_break_percentiles)),
                               labels = c("Low", "High")) +
            scale_y_continuous(breaks = as.numeric(quantile(p_leg$data[, 2],
                                                            legend_break_percentiles)),
                               labels = c("Low", "High")) +
            theme(axis.ticks = element_blank(),
                  axis.text.x = element_text(angle = 0, hjust = 0.5),
                  axis.title.x = element_text(vjust = 5),
                  axis.title.y = element_text(vjust = -5))
p_leg <- ggarrange(ggplot() + theme_void(), p_leg, ncol = 1)
p <- ggarrange(p_leg, p, widths = c(0.2, 0.8))

fig_name <- "commot_blended_plots"
file_name <- paste0(fig_paths[["main"]], fig_name,
                    output_figure_format)

ggsave(p, filename = file_name, units = "cm", dpi = 300, height = 15, width = 30)

12 Comparison of spatial characteristics

for (s in nice_samples) {
    c2l_cols <- grepl("c2l_", colnames(cells[[s]]@meta.data))
    c2l_md_mat <- as.matrix(cells[[s]]@meta.data[, c2l_cols])
    max_cts <- colnames(c2l_md_mat)[apply(c2l_md_mat, 1, which.max)]
    cells[[s]]$max_c2l_cts <- max_cts
}

plot_list <- list()
connections_df_list <- list()
for (s in nice_samples) {
    cells[[s]]$max_c2l_cts <- factor(cells[[s]]$max_c2l_cts, levels = colnames(c2l_md_mat))
    dists <- as.matrix(stats::dist(cells[[s]]@reductions$visual@cell.embeddings))
    imm_nbors_list <- lapply(1:ncol(cells[[s]]),
                             function(x) {
                                 imm_nbors_idxs <- dists[x, ] < 7.5 & dists[x, ] > 0
                                 imm_nbors_tab <- table(cells[[s]]$max_c2l_cts[imm_nbors_idxs])
                                 imm_nbors_tab <- as.numeric(imm_nbors_tab)

                                 return(imm_nbors_tab)
                             })
    imm_nbors_df <- do.call(rbind, imm_nbors_list)
    imm_nbors_df <- cbind(imm_nbors_df,
                          data.frame(cluster = cells[[s]]$max_c2l_cts))
    imm_nbors_df$cluster <- factor(imm_nbors_df$cluster,
                                   levels = levels(cells[[s]]$max_c2l_cts))

    connections_list <- lapply(levels(imm_nbors_df$cluster),
                               function(clus) {
                                   cluster_sub <- subset(imm_nbors_df, cluster == clus)
                                   if (nrow(cluster_sub) == 0) {
                                       cs <- rep(0, length(levels(imm_nbors_df$cluster)))
                                   } else {
                                       cs <- colSums(cluster_sub[-which(colnames(cluster_sub) == "cluster")])
                                   }
                                   names(cs) <- levels(imm_nbors_df$cluster)
                                   return(cs)
                               })
    names(connections_list) <- levels(imm_nbors_df$cluster)

    connections_mat_list <- lapply(names(connections_list),
                                   function(x) {
                                       mat <- as.matrix(connections_list[[x]])
                                       cbind(rownames(mat), mat, x)
                                   })
    connections_mat <- do.call(rbind, connections_mat_list)
    colnames(connections_mat) <- c("to", "num_links", "from")
    connections_df <- as.data.frame(connections_mat)
    connections_df$num_links <- as.numeric(connections_df$num_links)
    connections_df$from <- factor(connections_df$from,
                                  levels = levels(cells[[s]]$max_c2l_cts))
    connections_df$to <- factor(connections_df$to,
                                  levels = levels(cells[[s]]$max_c2l_cts))
    connections_df_list[[s]] <- connections_df
}

samples_connections_dims_list <- lapply(connections_df_list, function(x) {x$num_links})
samples_connections_dims <- do.call(rbind, samples_connections_dims_list)

samples_connections_dims <- t(apply(samples_connections_dims, 1, function(x) {x / sum(x)}))

samples_connections_dists <- as.matrix(stats::dist(samples_connections_dims))
samples_connections_dists <- max(samples_connections_dists) - samples_connections_dists
samples_connections_dists <- round(samples_connections_dists / max(samples_connections_dists), 2)
samples_connections_dists[upper.tri(samples_connections_dists)] <- NA
melted_samples_connections_dists <- melt(samples_connections_dists, na.rm = TRUE)
colnames(melted_samples_connections_dists)[3] <- "Similarity"

p <- ggplot(data = melted_samples_connections_dists,
            aes(x = Var1, y = Var2, fill = Similarity)) +
    geom_tile() +
    scale_x_discrete(labels = nice_samples) +
    scale_y_discrete(labels = nice_samples) +
    geom_text(aes(Var1, Var2, label = Similarity), color = "black", size = 4) +
    scale_fill_gradientn(colors = inferno(100)) +
    theme_minimal() +
    theme(axis.title = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1),
          axis.text = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_text(size = 14))

fig_name <- "spatial_connections_similarity_heatmap"
file_name <- paste0(fig_paths[["main"]], fig_name,
                    output_figure_format)
fig_list[["main"]][[fig_name]] <- p
ggsave(p, filename = file_name, units = "cm", height = 11.5, width = 15,
       dpi = 300)

samples_connections_dists_sub <- samples_connections_dists[4:6, 1:3]
estimated_pcws <- apply(samples_connections_dists_sub, 1,
                        function(x) {
                            c(x[1] * 7, x[2] * 11, x[3] * 17) / sum(x)
                        })
estimated_pcws <- data.frame(t(estimated_pcws))
estimated_pcws$org_name <- rownames(estimated_pcws)
estimated_pcws$org_name <- factor(estimated_pcws$org_name,
                                  levels = c("Donor1_day120_dif1_org1",
                                             "Donor2_day70_dif2_org2",
                                             "Donor1_day40_dif2_org2"))
estimated_pcws <- tidyr::pivot_longer(estimated_pcws,
                                      cols = -org_name,
                                      names_to = "tissue",
                                      values_to = "similarity")
estimated_pcws$tissue <- factor(estimated_pcws$tissue,
                                levels = c("PCW_17_REP_3", "PCW_11_REP_2",
                                           "PCW_7_REP_3"))

p <- ggplot(estimated_pcws,
            aes(estimated_pcws, fill = tissue, y = similarity,
                x = org_name)) +
        theme_bw() +
        geom_bar(position = "stack", stat = "identity") +
        scale_y_continuous(breaks = c(0, 7, 11, 17), limits = c(0, 17),
                           expand = expansion(mult = c(0, 0.02))) +
        scale_fill_manual(values=c("#D6618C", "#F9D15A", "#58A1E0"),
                          breaks = c(rev(levels(estimated_pcws$tissue)))) +
        geom_text(aes(label = after_stat(y), group = org_name),
                  stat = 'summary', fun = function(x) {round(sum(x), 1)},
                  hjust = -0.1) +
        coord_flip() +
        theme(axis.title.y = element_blank(),
              panel.grid.major.y = element_blank(),
              panel.grid.major.x = element_blank(),
              panel.grid.minor.x = element_blank(),
              legend.title = element_blank()) +
        ylab("Estimated PCW")
ggsave(p, filename = paste0("../paper_figs/main/estimated_pcw", output_figure_format),
       units = "cm", height = 4, width = 20, dpi = 300)

13 Assemble plots

# Figure 5
p <- wrap_plots(fig_list$fig_5[paste0("row_", 1:4)], ncol = 1,
                heights = c(3.5, 2, 2, 1.5))
file_name <- paste0(fig_paths[["main"]], "fig_5", output_figure_format)
ggsave(p, filename = file_name, units = "cm", height = 22,
       width = a4_width, dpi = 300)

# Figure 6
p <- ggarrange(plotlist = fig_list$fig_6[paste0("row_", 1:2)], ncol = 1,
               heights = c(0.3, 0.7))
file_name <- paste0(fig_paths[["main"]], "fig_6", output_figure_format)
# Need height of 18cm
ggsave(p, filename = file_name, units = "cm", height = 18,
       width = a4_width, dpi = 300)

14 Session info

Click to reveal details of R session

sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: aarch64-unknown-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/aarch64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/aarch64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] docknitr_1.0.1              rstudioapi_0.15.0          
##  [3] sys_3.4.2                   ggforce_0.4.2              
##  [5] showtext_0.9-7              showtextdb_3.0             
##  [7] sysfonts_0.8.9              msigdbr_7.5.1              
##  [9] knitr_1.45                  ggtext_0.1.2               
## [11] reshape2_1.4.4              tibble_3.2.1               
## [13] dplyr_1.1.3                 fgsea_1.24.0               
## [15] spacexr_2.2.1               scuttle_1.8.4              
## [17] SingleCellExperiment_1.20.1 DESeq2_1.38.3              
## [19] SummarizedExperiment_1.28.0 Biobase_2.58.0             
## [21] MatrixGenerics_1.10.0       matrixStats_1.0.0          
## [23] GenomicRanges_1.50.2        GenomeInfoDb_1.34.9        
## [25] IRanges_2.32.0              S4Vectors_0.36.2           
## [27] BiocGenerics_0.44.0         openxlsx_4.2.5.2           
## [29] ggrepel_0.9.4               ggpubr_0.6.0               
## [31] cowplot_1.1.1               ggplot2_3.4.4              
## [33] viridis_0.6.4               viridisLite_0.4.2          
## [35] patchwork_1.1.3             sctransform_0.4.1          
## [37] SeuratObject_5.0.0          Seurat_4.4.0               
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.4                spatstat.explore_3.2-5   
##   [3] reticulate_1.34.0         tidyselect_1.2.0         
##   [5] RSQLite_2.3.2             AnnotationDbi_1.60.2     
##   [7] htmlwidgets_1.6.2         BiocParallel_1.32.6      
##   [9] Rtsne_0.16                zellkonverter_1.8.0      
##  [11] munsell_0.5.0             codetools_0.2-19         
##  [13] ica_1.0-3                 future_1.33.0            
##  [15] miniUI_0.1.1.1            withr_2.5.2              
##  [17] spatstat.random_3.2-1     colorspace_2.1-0         
##  [19] progressr_0.14.0          filelock_1.0.2           
##  [21] highr_0.10                ROCR_1.0-11              
##  [23] ggsignif_0.6.4            tensor_1.5               
##  [25] listenv_0.9.0             labeling_0.4.3           
##  [27] GenomeInfoDbData_1.2.9    polyclip_1.10-6          
##  [29] farver_2.1.1              bit64_4.0.5              
##  [31] basilisk_1.10.2           parallelly_1.36.0        
##  [33] vctrs_0.6.4               generics_0.1.3           
##  [35] xfun_0.40                 Rnanoflann_0.0.2         
##  [37] markdown_1.12             doParallel_1.0.17        
##  [39] R6_2.5.1                  locfit_1.5-9.8           
##  [41] RcppZiggurat_0.1.6        hdf5r_1.3.10             
##  [43] bitops_1.0-7              spatstat.utils_3.0-4     
##  [45] cachem_1.0.8              DelayedArray_0.24.0      
##  [47] promises_1.2.1            scales_1.2.1             
##  [49] gtable_0.3.4              beachmat_2.14.2          
##  [51] globals_0.16.2            goftest_1.2-3            
##  [53] spam_2.10-0               rlang_1.1.1              
##  [55] splines_4.2.0             rstatix_0.7.2            
##  [57] lazyeval_0.2.2            spatstat.geom_3.2-7      
##  [59] broom_1.0.5               yaml_2.3.7               
##  [61] abind_1.4-5               backports_1.4.1          
##  [63] Rfast_2.1.0               httpuv_1.6.12            
##  [65] gridtext_0.1.5            tools_4.2.0              
##  [67] ellipsis_0.3.2            jquerylib_0.1.4          
##  [69] RColorBrewer_1.1-3        ggridges_0.5.4           
##  [71] Rcpp_1.0.11               plyr_1.8.9               
##  [73] sparseMatrixStats_1.10.0  zlibbioc_1.44.0          
##  [75] purrr_1.0.2               RCurl_1.98-1.12          
##  [77] basilisk.utils_1.10.0     deldir_1.0-9             
##  [79] pbapply_1.7-2             zoo_1.8-12               
##  [81] cluster_2.1.4             magrittr_2.0.3           
##  [83] data.table_1.14.8         scattermore_1.2          
##  [85] lmtest_0.9-40             RANN_2.6.1               
##  [87] fitdistrplus_1.1-11       mime_0.12                
##  [89] evaluate_0.22             xtable_1.8-4             
##  [91] XML_3.99-0.15             gridExtra_2.3            
##  [93] compiler_4.2.0            KernSmooth_2.23-22       
##  [95] crayon_1.5.2              htmltools_0.5.6.1        
##  [97] later_1.3.1               tidyr_1.3.0              
##  [99] geneplotter_1.76.0        Rfast2_0.1.5.2           
## [101] RcppParallel_5.1.7        DBI_1.1.3                
## [103] tweenr_2.0.3              MASS_7.3-60              
## [105] babelgene_22.9            Matrix_1.6-1.1           
## [107] car_3.1-2                 cli_3.6.1                
## [109] quadprog_1.5-8            parallel_4.2.0           
## [111] dotCall64_1.1-0           igraph_1.5.1             
## [113] pkgconfig_2.0.3           dir.expiry_1.6.0         
## [115] sp_2.1-1                  plotly_4.10.3            
## [117] spatstat.sparse_3.0-3     xml2_1.3.5               
## [119] foreach_1.5.2             annotate_1.76.0          
## [121] bslib_0.5.1               XVector_0.38.0           
## [123] stringr_1.5.0             digest_0.6.33            
## [125] RcppAnnoy_0.0.21          spatstat.data_3.0-3      
## [127] Biostrings_2.66.0         fastmatch_1.1-4          
## [129] rmarkdown_2.25            leiden_0.4.3             
## [131] uwot_0.1.16               DelayedMatrixStats_1.20.0
## [133] commonmark_1.9.0          shiny_1.7.5.1            
## [135] lifecycle_1.0.3           nlme_3.1-163             
## [137] jsonlite_1.8.7            carData_3.0-5            
## [139] fansi_1.0.5               pillar_1.9.0             
## [141] lattice_0.22-5            KEGGREST_1.38.0          
## [143] fastmap_1.1.1             httr_1.4.7               
## [145] survival_3.5-7            glue_1.6.2               
## [147] zip_2.3.0                 iterators_1.0.14         
## [149] png_0.1-8                 bit_4.0.5                
## [151] stringi_1.7.12            sass_0.4.7               
## [153] blob_1.2.4                memoise_2.0.1            
## [155] irlba_2.3.5.1             future.apply_1.11.0
installed.packages()[names(sessionInfo()$otherPkgs), "Version"]
##             docknitr           rstudioapi                  sys 
##              "1.0.1"             "0.15.0"              "3.4.2" 
##              ggforce             showtext           showtextdb 
##              "0.4.2"              "0.9-7"                "3.0" 
##             sysfonts              msigdbr                knitr 
##              "0.8.9"              "7.5.1"               "1.45" 
##               ggtext             reshape2               tibble 
##              "0.1.2"              "1.4.4"              "3.2.1" 
##                dplyr                fgsea              spacexr 
##              "1.1.3"             "1.24.0"              "2.2.1" 
##              scuttle SingleCellExperiment               DESeq2 
##              "1.8.4"             "1.20.1"             "1.38.3" 
## SummarizedExperiment              Biobase       MatrixGenerics 
##             "1.28.0"             "2.58.0"             "1.10.0" 
##          matrixStats        GenomicRanges         GenomeInfoDb 
##              "1.0.0"             "1.50.2"             "1.34.9" 
##              IRanges            S4Vectors         BiocGenerics 
##             "2.32.0"             "0.36.2"             "0.44.0" 
##             openxlsx              ggrepel               ggpubr 
##            "4.2.5.2"              "0.9.4"              "0.6.0" 
##              cowplot              ggplot2              viridis 
##              "1.1.1"              "3.4.4"              "0.6.4" 
##          viridisLite            patchwork          sctransform 
##              "0.4.2"              "1.1.3"              "0.4.1" 
##         SeuratObject               Seurat 
##              "5.0.0"              "4.4.0"